SU_25_River_Rockfish_Proj_Primary_V0.3

Setup, and Packages

Code Packages:

# Load packages
library(gridExtra) # For using grid.table to save tables as images

## Warning: package 'gridExtra' was built under R version 4.4.3

library(flextable) # Alternative method to save df as image

## Warning: package 'flextable' was built under R version 4.4.3

library(tidyverse)

## Warning: package 'ggplot2' was built under R version 4.4.3

## Warning: package 'lubridate' was built under R version 4.4.3

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.2     ✔ tibble    3.2.1
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::combine() masks gridExtra::combine()
## ✖ purrr::compose() masks flextable::compose()
## ✖ dplyr::filter()  masks stats::filter()
## ✖ dplyr::lag()     masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(zoo)

## 
## Attaching package: 'zoo'
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(lme4)

## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## 
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack

library(lmerTest)

## 
## Attaching package: 'lmerTest'
## 
## The following object is masked from 'package:lme4':
## 
##     lmer
## 
## The following object is masked from 'package:stats':
## 
##     step

library(rlang)

## 
## Attaching package: 'rlang'
## 
## The following objects are masked from 'package:purrr':
## 
##     %@%, flatten, flatten_chr, flatten_dbl, flatten_int, flatten_lgl,
##     flatten_raw, invoke, splice

library(nortest) # For ad.test
library(car) # For homogeneity of variance with leveneTest

## Warning: package 'car' was built under R version 4.4.3

## Loading required package: carData

## Warning: package 'carData' was built under R version 4.4.3

## 
## Attaching package: 'car'
## 
## The following object is masked from 'package:dplyr':
## 
##     recode
## 
## The following object is masked from 'package:purrr':
## 
##     some

library(rcompanion) # For scheirerRayHare test

## Warning: package 'rcompanion' was built under R version 4.4.3

library(pastecs) # For stat.desc

## Warning: package 'pastecs' was built under R version 4.4.3

## 
## Attaching package: 'pastecs'
## 
## The following objects are masked from 'package:dplyr':
## 
##     first, last
## 
## The following object is masked from 'package:tidyr':
## 
##     extract

library(ggpubr) # For assumption plots

## Warning: package 'ggpubr' was built under R version 4.4.3

## 
## Attaching package: 'ggpubr'
## 
## The following objects are masked from 'package:flextable':
## 
##     border, font, rotate

library(conover.test) # For the conover.Iman() test function
library(DescTools) # Contain alternative ConoverTest function

## Warning: package 'DescTools' was built under R version 4.4.3

## 
## Attaching package: 'DescTools'
## 
## The following object is masked from 'package:car':
## 
##     Recode

library(rstatix) # tidyverse adapted stat tests and mshapiro_test() functions

## Warning: package 'rstatix' was built under R version 4.4.3

## 
## Attaching package: 'rstatix'
## 
## The following object is masked from 'package:stats':
## 
##     filter

library(DHARMa) # For multi-level lm models and shapiro test

## Warning: package 'DHARMa' was built under R version 4.4.3

## This is DHARMa 0.4.7. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')

library(patchwork) # For combining ggplots

## Warning: package 'patchwork' was built under R version 4.4.3

Data Entry: Read-in data

# Read and call data into df

# Primary samples df
ABL90 <- read.csv("Raw_Data/ABL90_Raw.csv")

# River's Parturition metadata
partuition_subcat <- read.csv("Raw_Data/River's_Rockfish_Metadata_Parturition_V7.csv")

Data Filtering

# Data sifting: ABL90 dataset

# Step 1: Remove missing info
ABL_set1 <- ABL90 %>%
  filter(Patient.ID_edited != "") 
  

# Step 2: Separate Patient.ID into columns of sample types: blood plasma and instant freeze plasma 
ABL_set2 <- separate(ABL_set1, 'Patient.ID_edited', into = c("Patient.ID_edited", "Sample_Type"), sep = ",")  

## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 4 rows [489, 619, 1034,
## 1177].

# Step 3:
# Convert Patient.ID_edited column data to character type
ABL_set2$Patient.ID_edited <- as.character(ABL_set2$Patient.ID_edited) 

# Step 4 & 5: Sussing out specific sample errors  

# Remove insufficient samples
ABL_set3 <- ABL_set2 %>%
  filter(!is.na("Type")) %>%
  filter(!str_detect(Errors.detected.during.measurement, "Insufficient sample"))

# Remove inhomogeneous samples
ABL_set3 <- ABL_set3 %>%
  filter(!is.na("Type")) %>%
  filter(!str_detect(Errors.detected.during.measurement, "Inhomogeneous sample"))

# Step 6: Filter for only blood samples
ABL_b_samp <- ABL_set3 %>%
  filter(Sample_Type == "b")

Merge Datasets:

Join ABL90 df with parturition_subcat df: Making ABL_merged

# Rename columns to match: Change 'ID' in parturition_subcat to 'Patient.ID_edited' so it is ready to join with ABL90 df

# Please note the 'Treatment' col in parturition_subcat does not match with the finalized 'Ambient.Or.OAH' col in metadata_atresia_guide or ABL90_merged df

# Rename 'ID' col to 'patient.ID_edited'
partuition_subcat <- partuition_subcat %>%
  rename(Patient.ID_edited = ID)

# Connect main ABL90 dataset with my (River's) sub categorized parturition metadata 
# ABL_merged <- partuition_subcat %>%
# left_join(ABL_b_samp, by = "Patient.ID_edited")

ABL_merged <- ABL_b_samp %>%
inner_join(partuition_subcat, by = "Patient.ID_edited")

Rename Columns Neatly:

# Rename Treatment & Parturition Outcome
ABL_merged <- ABL_merged %>%
  rename(Ambient.Or.OAH = Consensus_General_Treatment,
         Pregnant.Or.Atresia = Consensus_Atresia_Or_Pregnant)

Remove Samples: Mortalities, NA’s, or Missing Info, and Replicates

Remove Replicate ID’s

# Cross Check all Columns: Cross validation to suss correct replicate row
check_params <- ABL_merged %>%
  select(Patient.ID_edited, Time, Sample.., Measuring.Mode, Ambient.Or.OAH, Pregnant.Or.Atresia, pH, Glu..mmol.L., Hct...., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L., pCO2..mmHg.) %>%
  filter(Patient.ID_edited == "9782D")
print(check_params)

##   Patient.ID_edited             Time Sample.. Measuring.Mode Ambient.Or.OAH
## 1             9782D 12/28/2023 10:04     1016         C 65uL        Ambient
## 2             9782D   1/9/2024 17:21      863         S 65uL        Ambient
##   Pregnant.Or.Atresia    pH Glu..mmol.L. Hct.... Na...mmol.L. Cl...mmol.L.
## 1    Post Parturition 7.348          0.5     9.8          184          149
## 2    Post Parturition 7.446          0.7    13.3          181          161
##   K...mmol.L. Ca....mmol.L. pCO2..mmHg.
## 1         2.8          1.32         4.2
## 2         2.9          1.31         6.8

check_params <- ABL_merged %>%
  select(Patient.ID_edited, Time, Sample.., Measuring.Mode, Ambient.Or.OAH, Pregnant.Or.Atresia, pH, Glu..mmol.L., Hct...., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L., pCO2..mmHg.) %>%
  filter(Patient.ID_edited == "9783D")
print(check_params)

##   Patient.ID_edited            Time Sample.. Measuring.Mode Ambient.Or.OAH
## 1             9783D 1/19/2024 10:37      934         S 65uL        Ambient
## 2             9783D  1/7/2024 12:32      860         S 65uL        Ambient
##   Pregnant.Or.Atresia    pH Glu..mmol.L. Hct.... Na...mmol.L. Cl...mmol.L.
## 1    Post Parturition 7.484          1.7     0.2          175          193
## 2    Post Parturition 7.450          1.3    14.7          179          164
##   K...mmol.L. Ca....mmol.L. pCO2..mmHg.
## 1        34.5          1.29         8.5
## 2         2.5          1.30         5.3

# Removed replicates: 9782D (2x) and 9783D (2x)
ABL_merged <- ABL_merged %>%
  filter(Sample.. != "863",
         Sample.. != "934")

Review Replicates:

Row	ID	Time	Sample..	Measuring.Mode	pH	Glu.mmol.L.	Hct….	Na…mmol.L.	Cl…mmol.L.	K…mmol.L.	Ca.mmol.L.
32	9782D	12/28/2023 10:04	1016	C 65uL	7.348	0.5	9.8	184	149	2.8	1.32
39	9782D	1/9/2024 17:21	863	S 65uL	7.446	0.7	13.3	181	161	2.9	1.31
36	9783D	1/19/2024 10:37	934	S 65uL	7.484	1.7	0.2	175	193	34.5	1.29
41	9783D	1/7/2024 12:32	860	S 65uL	7.450	1.3	14.7	179	164	2.5	1.30

Methods for replicate sample removal: Check processing date of IDs, and remove samples that do not match that date.

Case ID_9782D:

• Origianl ABL90 set contains:

• 9782D,b (Sample.. 1016, Time = 12/28/2023 10:04:00 AM)

• 9782D,p (Sample.. 1020, Time = 12/28/2023 10:25:00 AM)

• 9782D,p (Sample.. 865, Time = 1/9/2024 5:29:00 PM)

• 9782D,b (Sample.. 863, Time = 1/9/2024 5:21:00 PM)

• 2024_Fish_Metadata.csv says ID_9782D was processed 3/29/2024

• Note Processing date in metadata will not match time date stamp in ABL90 data due to the machine being reset, and innacurrate dates being automatically recorded.

Deciding info: - For ID-9782D, refer to metadata sheet and look for other samples with that same processing date, then view thoes samples ABL90 ‘Time’. Try to find a matching ‘Time’ with other IDs and 9782D that were processed together during the same time recorded by the machine.

• ID_97838,b (single obs, Sample.. 1024) Time = 12/28/2023 11:49:00 AM

• Since 9782D,b Sample.. 1016 and 97838,b Sample.. 1024 both have a date timestamp occurrence on 12/28/2023, I have decided to keep 9782D,b Sample.. 1016 and drop 9783D,b Sample.. 863

Case ID_9783D

• Original ABL90 dataset contains:

• 9783D,b (Sample.. 934, Time = 1/19/2024 10:37:00 AM)

• 9783D,c (Sample.. 861, Time = 1/7/2024 12:38:00 PM)

• 9783D,b (Sample.. 860, Time = 1/7/2024 12:32:00 PM)

• 2024_Fish_Metadata.csv says ID_9783D was processed 3/20/2024

• Informed that ID_9783D: Sample.. 934 is a plasma sample (on account of its low hct value), remove this and keep Sample.. 860 (parameter values appear correct).

Samples: Data subsets

Remove Motalities and No info IDs and assign analysis ready data-frames:

# New df with Moralities removed: Note none of these samples made it into ABL90 df anyway, so looks like they are already filtered out.
 Live_Samples <- ABL_merged %>%
  filter(Patient.ID_edited != "9780C", # Mortality
         Patient.ID_edited != "777AE", # Mortality
         Patient.ID_edited != "777CA") # Mortality after parturition


# New df with mortality and 'No info' Id's removed
 Primary_Samples <- Live_Samples %>%
   filter(Patient.ID_edited != "777A0", # No info
           Patient.ID_edited != "9782F", # No info
           Patient.ID_edited != "777B3", # No info
           Patient.ID_edited != "777AA", # No info
           Patient.ID_edited != "777DE", # No info
           Patient.ID_edited != "777CE", # No info
           Patient.ID_edited != "777A6") # No info
 
 
# New df of Only Ambient Treatment: For testing parturition success
 Ambient_Samples <- Primary_Samples %>%
   filter(Ambient.Or.OAH == "Ambient")
 
# New df of Fecundity Samples:
 Fecundity_Samples <- Primary_Samples %>%
   filter(Fecundity_Or_Brood_Count != "NA",
     Fecundity_Class != "NA") # removes 97706

# Fecundity samples without atresia Ids
 Fecundity_No_Atresia_Samples <- Primary_Samples %>%
   filter(All_Fecundity != "Na",
          All_Fecundity != 0)  # removes atresia IDs

# Fecundity Ambient Samples
 Amb_Fec_Samples <- Primary_Samples %>%
   filter(Ambient.Or.OAH == "Ambient",
          All_Fecundity != "NA") # removes OAH treatment samples and ID 97706
 
# Fecundity Ambient Samples without atresia Ids
 Amb_Fec_No_Atresia_Samples <- Primary_Samples %>%
   filter(Ambient.Or.OAH == "Ambient", # removes OAH treatment 
          All_Fecundity != "NA", # removes all missing info IDs
          All_Fecundity != 0) # removes atresia IDs

Change Data Types:

# Change Continuous Data to type to numeric, double, or factor

# Change data type of ions to factor or numeric
#glimpse(General_Samples)
#glimpse(Ambient_Samples)
 
Primary_Samples <- Primary_Samples %>%
  mutate(across(c(All_Fecundity, pH, Hct...., Glu..mmol.L., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L.), as.numeric))

Ambient_Samples <- Ambient_Samples %>%
  mutate(across(c(All_Fecundity, pH, Hct...., Glu..mmol.L., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L.), as.numeric))

Fecundity_Samples <- Fecundity_Samples %>%
  mutate(across(c(All_Fecundity, pH, Hct...., Glu..mmol.L., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L.), as.numeric))

Fecundity_No_Atresia_Samples <- Fecundity_Samples %>%
  mutate(across(c(All_Fecundity, pH, Hct...., Glu..mmol.L., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L.), as.numeric))

Amb_Fec_Samples <- Amb_Fec_Samples %>%
  mutate(across(c(All_Fecundity, pH, Hct...., Glu..mmol.L., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L.), as.numeric))
  
Amb_Fec_No_Atresia_Samples <- Amb_Fec_No_Atresia_Samples %>%
  mutate(across(c(All_Fecundity, pH, Hct...., Glu..mmol.L., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L.), as.numeric))

# glimpse(Primary_Samples)
# glimpse(Ambient_Samples)
# glimpse(Fecundity_Samples)
# glimpse(Amb_Fec_Samples)

Order Factor Levels

# Fecundity data 

#unique(Fecundity_Samples$Fecundity_Class)

# Arrange the order of parturition categories for visualizations & tests

# Primary Samples

# Arrange Treatment
Primary_Samples$Ambient.Or.OAH <- ordered(Primary_Samples$Ambient.Or.OAH, levels = c("Ambient", "OAH"))

# Arrange Parturition outcome
Primary_Samples$Pregnant.Or.Atresia <- ordered(Primary_Samples$Pregnant.Or.Atresia, levels = c("Post Parturition", "Atresia"))

# Arrange Brood Condition
Primary_Samples$Consensus_Brood_Condition <- ordered(Primary_Samples$Consensus_Brood_Condition, levels = c("Excellent", "Normal", "Low", "Very Low", "Atresia"))

# Arrange Fecundity Class
Primary_Samples$Fecundity_Class <- ordered(Primary_Samples$Fecundity_Class, levels = c("High (>50,000)", "Low (~1,000s)", "Atresia"))

# Remove NAs from ordered character factors
Primary_Samples <- Primary_Samples %>%
  filter(Fecundity_Class != is.na(Fecundity_Class)) %>%
  filter(Consensus_Brood_Condition != is.na(Consensus_Brood_Condition))
  

#class(Primary_Samples$Fecundity_Class)

Save filtered dataset into data-worked folder

# Save merged dataset in a worked folder

write.csv(x = Primary_Samples, file = "Worked-Data/Primary_Samples")

write.csv(x = Ambient_Samples, file = "Worked-Data/Ambient_Samples")

write.csv(x = Fecundity_Samples, file = "Worked-Data/Fecundity_Samples")

write.csv(x = Fecundity_No_Atresia_Samples, file = "Worked-Data/Fecundity_No_Atresia_Samples")

write.csv(x = Amb_Fec_Samples, file = "Worked-Data/Amb-Fec_Samples")

write.csv(x = Amb_Fec_No_Atresia_Samples, file = "Worked-Data/Amb_Fec_No_Atresia_Samples")

Sample Size: n

# Sample Size

# Primary data: Treatment & Pregnancy outcome
primary_sample_table <- Primary_Samples %>%
  group_by(Ambient.Or.OAH, Pregnant.Or.Atresia) %>%
  summarise(n_size = n_distinct(Patient.ID_edited), .groups = "drop")
print(primary_sample_table)

## # A tibble: 4 × 3
##   Ambient.Or.OAH Pregnant.Or.Atresia n_size
##   <ord>          <ord>                <int>
## 1 Ambient        Post Parturition         8
## 2 Ambient        Atresia                 13
## 3 OAH            Post Parturition         2
## 4 OAH            Atresia                  5

 pdf("data-images/primary_sample_table.pdf")
 grid.table(primary_sample_table)
 dev.off()

## png 
##   2

 png("data-images/primary_sample_table.png")
 grid.table(primary_sample_table)
 dev.off()

## png 
##   2

# Primary data: Brood Condition
primary_sample_table_brood_codnition <- Primary_Samples %>%
  group_by(Ambient.Or.OAH, Consensus_Brood_Condition) %>%
  summarise(n_size = n_distinct(Patient.ID_edited), .groups = "drop")
print(primary_sample_table_brood_codnition)

## # A tibble: 6 × 3
##   Ambient.Or.OAH Consensus_Brood_Condition n_size
##   <ord>          <ord>                      <int>
## 1 Ambient        Excellent                      1
## 2 Ambient        Normal                         4
## 3 Ambient        Low                            3
## 4 Ambient        Atresia                       13
## 5 OAH            Excellent                      2
## 6 OAH            Atresia                        5

# Primary data: Fecundity Class
primary_sample_table_fecundity_class <- Primary_Samples %>%
  group_by(Ambient.Or.OAH, Pregnant.Or.Atresia, Fecundity_Class) %>%
  summarise(n_size = n_distinct(Patient.ID_edited), .groups = "drop")
print(primary_sample_table_fecundity_class)

## # A tibble: 5 × 4
##   Ambient.Or.OAH Pregnant.Or.Atresia Fecundity_Class n_size
##   <ord>          <ord>               <ord>            <int>
## 1 Ambient        Post Parturition    High (>50,000)       4
## 2 Ambient        Post Parturition    Low (~1,000s)        4
## 3 Ambient        Atresia             Atresia             13
## 4 OAH            Post Parturition    High (>50,000)       2
## 5 OAH            Atresia             Atresia              5

Sample Size Barplot:

# View Sample Distributions

# Sample Size barplot

# Primary Samples: General sample size 
primary_sample_barplot <- Primary_Samples %>%
  group_by(Pregnant.Or.Atresia, Ambient.Or.OAH) %>% 
  summarise(n_size = n_distinct(Patient.ID_edited), .groups = "drop") %>%
  ggplot(aes(x = Pregnant.Or.Atresia, y = n_size)) +
  geom_bar(stat = "identity" ) +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 15.5), breaks = c(0, 3, 6, 9, 12, 15)) +
  facet_grid(~ Ambient.Or.OAH) +
  labs(title = "Primary Sample Size",
        x = "Parturition Status & Treatment Type",
        y = "Sample Size (n)") +
  guides(fill = guide_legend((title = "Reproductive State"))) +
  theme_classic() +
  theme(plot.title = element_text(size = 30, hjust = 0.5))

print(primary_sample_barplot)

ggsave(filename = "data-images/.png", plot = primary_sample_barplot, device = "png")

## Saving 7 x 5 in image

ggsave(filename = "data-images/primary_sample_barplot.pdf", plot = primary_sample_barplot, device = "pdf")

## Saving 7 x 5 in image

# Primary samples: Fecundity Class sample size
primary_sample_fecundity_class_barplot <- Primary_Samples %>%
  group_by(Fecundity_Class, Ambient.Or.OAH) %>% 
  summarise(n_size = n_distinct(Patient.ID_edited), .groups = "drop") %>%
  ggplot(aes(x = Fecundity_Class, y = n_size)) +
  geom_bar(stat = "identity" ) +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 15.5), breaks = c(0, 3, 6, 9, 12, 15)) +
  facet_grid(~ Ambient.Or.OAH) +
  labs(title = "Sample Size",
        x = "Parturition Status & Treatment Type",
        y = "Sample Size (n)") +
  guides(fill = guide_legend((title = "Reproductive State"))) +
  theme_classic() +
  theme(plot.title = element_text(size = 30, hjust = 0.5))

print(primary_sample_fecundity_class_barplot)

# Primary samples: Brood Condition sample size
primary_sample_brood_condition_barplot <- Primary_Samples %>%
  group_by(Consensus_Brood_Condition, Ambient.Or.OAH) %>% 
  summarise(n_size = n_distinct(Patient.ID_edited), .groups = "drop") %>%
  ggplot(aes(x = Consensus_Brood_Condition, y = n_size)) +
  geom_bar(stat = "identity" ) +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 15.5), breaks = c(0, 3, 6, 9, 12, 15)) +
  facet_grid(~ Ambient.Or.OAH) +
  labs(title = "Sample Size",
        x = "Parturition Status & Treatment Type",
        y = "Sample Size (n)") +
  guides(fill = guide_legend((title = "Reproductive State"))) +
  theme_classic() +
  theme(plot.title = element_text(size = 30, hjust = 0.5))

print(primary_sample_brood_condition_barplot)

pH

pH Summary Stats

# pH
# Summary stats

# Primary data
Primary_Samples %>%
  group_by(Pregnant.Or.Atresia, Ambient.Or.OAH) %>%
  summarize(count = n(),
            median = round(median(pH), 3),
            mean = round(mean(pH), 3),
            sd = round(sd(pH), 3),
            cv = round(sd(pH)/mean(pH), 3)) %>%
  ungroup()

## `summarise()` has grouped output by 'Pregnant.Or.Atresia'. You can override
## using the `.groups` argument.

## # A tibble: 4 × 7
##   Pregnant.Or.Atresia Ambient.Or.OAH count median  mean    sd    cv
##   <ord>               <ord>          <int>  <dbl> <dbl> <dbl> <dbl>
## 1 Post Parturition    Ambient            8   7.44  7.44 0.063 0.008
## 2 Post Parturition    OAH                2   7.45  7.45 0.057 0.008
## 3 Atresia             Ambient           13   7.37  7.38 0.1   0.014
## 4 Atresia             OAH                5   7.24  7.25 0.058 0.008

pH Boxplot

# pH boxplot

# Primary samples
pH_primary_boxplot <- ggplot(data = Primary_Samples, aes(x = Pregnant.Or.Atresia, y = pH, fill = Ambient.Or.OAH)) +
  geom_boxplot(aes(x = Pregnant.Or.Atresia, y = pH, fill = Ambient.Or.OAH)) +
  geom_point(aes(color = Ambient.Or.OAH), position = position_dodge(width = 0.75), alpha = 0.3, color = "black") +
  labs(title = "Primary pH", x = "Parturition Type", y = "pH") +
  guides(fill = guide_legend((title = "Treatment Type"))) + 
  scale_fill_manual(values = c("Ambient" = "#189bff", "OAH" = "#598c78")) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(pH_primary_boxplot)

ggsave(filename = "data-images/pH_primary_boxplot.pdf", plot = pH_primary_boxplot, device = "pdf")

## Saving 7 x 5 in image

ggsave(filename = "data-images/pH_primary_boxplot.png", plot = pH_primary_boxplot, device = "png")

## Saving 7 x 5 in image

pH Scatterplot: No Atresia IDs

# pH
# Primary Samples

# Weight adjusted scatterplot 

# Primary scatterplot: No atresia ids 
pH_primary_weight_adj_scatterplot <- ggplot(data = Primary_Samples) + 
  geom_point(aes(x = pH, y = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight), colour = "black") +
  labs(title = "pH", x = "pH", y = "Weight Adjusted Fecundity") +
  scale_y_continuous() + 
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(pH_primary_weight_adj_scatterplot)

## Warning: Removed 18 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggsave(filename = "data-images/pH_primary_weight_adj_scatterplot.pdf", plot = pH_primary_weight_adj_scatterplot, device = "pdf")

## Saving 7 x 5 in image

## Warning: Removed 18 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggsave(filename = "data-images/pH_primary_weight_adj_scatterplot.png", plot = pH_primary_weight_adj_scatterplot, device = "png")

## Saving 7 x 5 in image

## Warning: Removed 18 rows containing missing values or values outside the scale range
## (`geom_point()`).

# Facet by Treatment
pH_primary_weight_adj_scatterplot <- ggplot(data = Primary_Samples) + 
  geom_point(aes(x = pH, y = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight), colour = "black") +
  labs(title = "pH", x = "pH", y = "Weight Adjusted Fecundity") +
  scale_y_continuous() + 
  theme_classic() +
  facet_wrap(~ Ambient.Or.OAH) +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(pH_primary_weight_adj_scatterplot)

## Warning: Removed 18 rows containing missing values or values outside the scale range
## (`geom_point()`).

Data Distribution

# pH
# Parametric Assumptions

# Data Distribution
hist(Primary_Samples$pH)

plotNormalHistogram(Primary_Samples$pH) 

plotNormalDensity(Primary_Samples$pH)

pH Stat Tests

# pH
# Primary Samples

# Interactive aov model
pH_primary_aov_int <- aov(pH ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)
summary(pH_primary_aov_int) # interaction almost significant

##                                    Df  Sum Sq Mean Sq F value  Pr(>F)   
## Pregnant.Or.Atresia                 1 0.06629 0.06629   9.638 0.00483 **
## Ambient.Or.OAH                      1 0.03877 0.03877   5.637 0.02593 * 
## Pregnant.Or.Atresia:Ambient.Or.OAH  1 0.02078 0.02078   3.022 0.09495 . 
## Residuals                          24 0.16506 0.00688                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Additive aov model
# If interaction not significant (p-value > 0.05), use additive model
# pH_primary_aov_additive <- aov(pH ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = Primary_Samples)
# summary(pH_primary_aov_additive)

# Parametric Post-Hoc: Tukey's HSD test
TukeyHSD(pH_primary_aov_int)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = pH ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)
## 
## $Pregnant.Or.Atresia
##                                diff        lwr         upr    p adj
## Atresia-Post Parturition -0.1015444 -0.1690501 -0.03403877 0.004834
## 
## $Ambient.Or.OAH
##                    diff        lwr         upr    p adj
## OAH-Ambient -0.08561481 -0.1603143 -0.01091529 0.026428
## 
## $`Pregnant.Or.Atresia:Ambient.Or.OAH`
##                                                      diff        lwr
## Atresia:Ambient-Post Parturition:Ambient      -0.06418269 -0.1669825
## Post Parturition:OAH-Post Parturition:Ambient  0.00862500 -0.1722337
## Atresia:OAH-Post Parturition:Ambient          -0.19247500 -0.3228940
## Post Parturition:OAH-Atresia:Ambient           0.07280769 -0.1009557
## Atresia:OAH-Atresia:Ambient                   -0.12829231 -0.2486791
## Atresia:OAH-Post Parturition:OAH              -0.20110000 -0.3925028
##                                                       upr     p adj
## Atresia:Ambient-Post Parturition:Ambient       0.03861711 0.3345551
## Post Parturition:OAH-Post Parturition:Ambient  0.18948366 0.9991631
## Atresia:OAH-Post Parturition:Ambient          -0.06205597 0.0023249
## Post Parturition:OAH-Atresia:Ambient           0.24657107 0.6594942
## Atresia:OAH-Atresia:Ambient                   -0.00790551 0.0337411
## Atresia:OAH-Post Parturition:OAH              -0.00969719 0.0369572

# Normality plot

# Qqplot: Residuals of aov or lm model (plot residuals of aov vs lm are identical)
pH_primary_res_qqplot <- ggqqplot(residuals(aov(pH ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples))) +
labs(title = "Primary pH Interactive Residual QQplot",
     subtitle = "residuals(aov(pH ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples))",
                               x = "pH Theoretical Quantiles (Predicted values)",
                               y = "pH Sample Quantiles") +
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5))
print(pH_primary_res_qqplot)

# Normality Test
shapiro.test(Primary_Samples$pH)

## 
##  Shapiro-Wilk normality test
## 
## data:  Primary_Samples$pH
## W = 0.97662, p-value = 0.7632

shapiro.test(residuals(lm(pH ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(lm(pH ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples))
## W = 0.96397, p-value = 0.431

# Parametric variance test:
bartlett.test(pH ~ Pregnant.Or.Atresia, data = Primary_Samples) # nonparametric

## 
##  Bartlett test of homogeneity of variances
## 
## data:  pH by Pregnant.Or.Atresia
## Bartlett's K-squared = 3.4356, df = 1, p-value = 0.0638

bartlett.test(pH ~ Ambient.Or.OAH, data = Primary_Samples)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  pH by Ambient.Or.OAH
## Bartlett's K-squared = 0.33993, df = 1, p-value = 0.5599

# Nonparametric variance test:
LeveneTest(pH ~ Pregnant.Or.Atresia, data = Primary_Samples) 

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  2.3324 0.1388
##       26

leveneTest(pH ~ Ambient.Or.OAH, data = Primary_Samples)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.1912 0.6655
##       26

# Nonparametric Stat test & Post-Hoc:

kruskal.test(pH ~ Pregnant.Or.Atresia, data = Primary_Samples) # Significant

## 
##  Kruskal-Wallis rank sum test
## 
## data:  pH by Pregnant.Or.Atresia
## Kruskal-Wallis chi-squared = 6.954, df = 1, p-value = 0.008363

kruskal.test(pH ~ Ambient.Or.OAH, data = Primary_Samples) # Significant

## 
##  Kruskal-Wallis rank sum test
## 
## data:  pH by Ambient.Or.OAH
## Kruskal-Wallis chi-squared = 3.7502, df = 1, p-value = 0.0528

DunnTest(pH ~ Pregnant.Or.Atresia, data = Primary_Samples)

## 
##  Dunn's test of multiple comparisons using rank sums : holm  
## 
##                          mean.rank.diff   pval    
## Atresia-Post Parturition      -8.555556 0.0084 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

DunnTest(pH ~ Ambient.Or.OAH, data = Primary_Samples)

## 
##  Dunn's test of multiple comparisons using rank sums : holm  
## 
##             mean.rank.diff   pval    
## OAH-Ambient      -6.952381 0.0528 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Parametric Results: Ambient Fecundity No Atresia Sample

• A significant difference detected in pH between pregnancy groups and treatment groups

• Interaction almost significant (), interactive model used

• Parametric post Hoc: Tukey HSD with interactive model

• Ask how to interpret interactive results…

Parametric Assumptions: Across Sample comparisons…

• Sample size (n > 30): Not Met

• Normality: Met

• Homoscedasticity: Not Met, nonparametric

• Evidence: data: pH by Pregnant.Or.Atresia Bartlett’s K-squared = 3.4356, df = 1, p-value = 0.0638

• Nonparametric variance test (Levene’s test): Met

Non-Parametric Results: - A significant increase detected in pH for post parturition group compared to atresia group. - Evidence: Dunn’s test pval = 0.0084 **

• A significant increase detected in pH for ambient group compared to OAH group.
• Evidence: Dunn’s test pval = 0.0528

# ph
# Primary Samples

# Regression Model 

# Install flexplot
# devtools::install_github("dustinfife/flexplot", ref="development")
require(flexplot)

## Loading required package: flexplot

## 
## Attaching package: 'flexplot'

## The following object is masked from 'package:ggplot2':
## 
##     flip_data

# If interaction not significant, use additive model
# p-value significant if (a < 0.1 or a < 0.05)
pH_primary_lm <- lm(pH ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)
summary(pH_primary_lm)

## 
## Call:
## lm(formula = pH ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.149692 -0.042344 -0.008546  0.046671  0.222308 
## 
## Coefficients:
##                                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                             7.37887    0.01969 374.762  < 2e-16 ***
## Pregnant.Or.Atresia.L                  -0.09379    0.02785  -3.368  0.00255 ** 
## Ambient.Or.OAH.L                       -0.04231    0.02785  -1.519  0.14172    
## Pregnant.Or.Atresia.L:Ambient.Or.OAH.L -0.06846    0.03938  -1.738  0.09495 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.08293 on 24 degrees of freedom
## Multiple R-squared:  0.4326, Adjusted R-squared:  0.3617 
## F-statistic: 6.099 on 3 and 24 DF,  p-value: 0.003094

pH_primary_lm_res_plot <- simulateResiduals(fittedModel = pH_primary_lm)
plot(pH_primary_lm_res_plot)

visualize(pH_primary_lm, plot = "model")

# Residuals of model for shapiro test
pH_primary_lm_res <- residuals(pH_primary_lm)

# Normality test
shapiro.test(pH_primary_lm_res)

## 
##  Shapiro-Wilk normality test
## 
## data:  pH_primary_lm_res
## W = 0.96397, p-value = 0.431

# Scedasticity test: 
# bartlett's test parametric
bartlett.test(pH ~ Pregnant.Or.Atresia, data = Primary_Samples) # nonparametric

## 
##  Bartlett test of homogeneity of variances
## 
## data:  pH by Pregnant.Or.Atresia
## Bartlett's K-squared = 3.4356, df = 1, p-value = 0.0638

bartlett.test(pH ~ Ambient.Or.OAH, data = Primary_Samples)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  pH by Ambient.Or.OAH
## Bartlett's K-squared = 0.33993, df = 1, p-value = 0.5599

# Levene's test if nonparametric
leveneTest(pH ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3  0.6518 0.5896
##       24

leveneTest(pH ~ Pregnant.Or.Atresia, data = Primary_Samples) # Close to significant

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  2.3324 0.1388
##       26

leveneTest(pH ~ Ambient.Or.OAH, data = Primary_Samples)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.1912 0.6655
##       26

# Non-parametric regression test???:

Results: Primary Sample (With both treatments and parturition groups)

Linear model interaction significant (Pr(>|t|) = 0.09495 .), interactive model used.

Linear Model: Interactive lm param ~ parturition * treatment - A significant relationship detected between pH, parturition groups, and treatment - Evidence: summary lm Call: lm(formula = pH ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples) Residual standard error: 0.08293 on 24 degrees of freedom Multiple R-squared: 0.4326, Adjusted R-squared: 0.3617 F-statistic: 6.099 on 3 and 24 DF, p-value: 0.003094

• How to know which groups are correlated, and how strongly?
• (test lm individually)

Parametric Assumptions: Across Sample comparisons…

• Sample size (n > 30): Not Met

• Normality (of residuals) : Met

• Homoscedasticity: Met

# ph
# Primary Samples

# Regression Model pt 2

# Individually assess parameter to parturition groups & treatment groups
pH_primary_lm1 <- lm(pH ~ Pregnant.Or.Atresia, data = Primary_Samples)

pH_primary_lm2 <- lm(pH ~ Ambient.Or.OAH, data = Primary_Samples)

# p-value significant if (a < 0.1 or a < 0.05)
summary(pH_primary_lm1) # Significant

## 
## Call:
## lm(formula = pH ~ Pregnant.Or.Atresia, data = Primary_Samples)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.152056 -0.057714  0.000944  0.042900  0.257944 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            7.39183    0.01833  403.29   <2e-16 ***
## Pregnant.Or.Atresia.L -0.07180    0.02592   -2.77   0.0102 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09295 on 26 degrees of freedom
## Multiple R-squared:  0.2279, Adjusted R-squared:  0.1982 
## F-statistic: 7.673 on 1 and 26 DF,  p-value: 0.01021

summary(pH_primary_lm2) # Significant

## 
## Call:
## lm(formula = pH ~ Ambient.Or.OAH, data = Primary_Samples)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.17414 -0.06893 -0.01614  0.06336  0.19786 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       7.35350    0.02111 348.403   <2e-16 ***
## Ambient.Or.OAH.L -0.06738    0.02985  -2.257   0.0326 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09672 on 26 degrees of freedom
## Multiple R-squared:  0.1639, Adjusted R-squared:  0.1317 
## F-statistic: 5.095 on 1 and 26 DF,  p-value: 0.03262

pH_primary_lm_res_plot1 <- simulateResiduals(fittedModel = pH_primary_lm1)
plot(pH_primary_lm_res_plot1)

pH_primary_lm_res_plot2 <- simulateResiduals(fittedModel = pH_primary_lm2)
plot(pH_primary_lm_res_plot2)

# For shapiro test
pH_primary_lm_res1 <- residuals(pH_primary_lm1)

pH_primary_lm_res2 <- residuals(pH_primary_lm2)

# Normality test
shapiro.test(pH_primary_lm_res1) # Normal

## 
##  Shapiro-Wilk normality test
## 
## data:  pH_primary_lm_res1
## W = 0.96109, p-value = 0.3699

shapiro.test(pH_primary_lm_res2) # Normal

## 
##  Shapiro-Wilk normality test
## 
## data:  pH_primary_lm_res2
## W = 0.97843, p-value = 0.8112

Results: Primary Sample (With both treatments and parturition groups)

Linear Model: Individual test correlation between param and treatment, and param and parturition outcome

Correlation: pH and Post-parturition vs atresia - A significant relationship detected with pH and parturition groups - Evidence: Residual standard error: 0.09295 on 26 degrees of freedom Multiple R-squared: 0.2279, Adjusted R-squared: 0.1982 F-statistic: 7.673 on 1 and 26 DF, p-value: 0.01021

• Meaning: Post parturition groups display a significant correlative increase in pH, compared to the atresia group (p-value: 0.01021).

Warning: DHARMA states no significant (n.s.) within group deviation from uniformity, and the Levene Test for homogeneity of variance n.s.

Parametric Assumptions: Across Sample comparisons…

• Sample size (n > 30): Not Met

• Normality: Met

• Homoscedasticity: Met

Residuals Plot (if significant)

# pH
# Primary samples

# ggplot with geom_smooth(aes(x=,y=)), method = "lm")


# lm for pH & parturition scatterplot: 
pH_primary_lm_scatterplot1 <- ggplot(data = Primary_Samples, aes(x = Pregnant.Or.Atresia, y = pH)) + 
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "pH", x = "Parturition Group", y = "pH", color = "Category") +
  scale_y_continuous() + 
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(pH_primary_lm_scatterplot1)

## `geom_smooth()` using formula = 'y ~ x'

ggsave(filename = "data-images/pH_primary_lm_scatterplot1.pdf", plot = pH_primary_lm_scatterplot1, device = "pdf")

## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'

ggsave(filename = "data-images/pH_primary_lm_scatterplot1.png", plot = pH_primary_lm_scatterplot1, device = "png")

## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'

# lm for pH & parturition scatterplot: faceted by treatment
pH_primary_lm_scatterplot1.1 <- ggplot(data = Primary_Samples, aes(x = Pregnant.Or.Atresia, y = pH)) + 
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "pH", x = "Parturition Group", y = "pH", color = "Category") +
  scale_y_continuous() + 
  theme_classic() +
  facet_wrap(~ Ambient.Or.OAH) +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(pH_primary_lm_scatterplot1.1)

## `geom_smooth()` using formula = 'y ~ x'

# lm for pH & treatment scatterplot:
pH_primary_lm_scatterplot2 <- ggplot(data = Primary_Samples, aes(x = Ambient.Or.OAH, y = pH)) + 
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "pH", x = "Treatment Group", y = "pH", color = "Category") +
  scale_y_continuous() + 
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(pH_primary_lm_scatterplot2)

## `geom_smooth()` using formula = 'y ~ x'

ggsave(filename = "data-images/pH_primary_lm_scatterplot2.pdf", plot = pH_primary_lm_scatterplot2, device = "pdf")

## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'

ggsave(filename = "data-images/pH_primary_lm_scatterplot2.png", plot = pH_primary_lm_scatterplot2, device = "png")

## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'

# lm for pH & treatment scatterplot: faceted by parturition
pH_primary_lm_scatterplot2.1 <- ggplot(data = Primary_Samples, aes(x = Ambient.Or.OAH, y = pH)) + 
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "pH", x = "Treatment Group", y = "pH", color = "Category") +
  scale_y_continuous() + 
  theme_classic() +
  facet_wrap(~ Pregnant.Or.Atresia) +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(pH_primary_lm_scatterplot2.1)

## `geom_smooth()` using formula = 'y ~ x'

# pH
# lm Boxplot

pH_primary_lm_boxplot1 <- ggplot(data = Primary_Samples, aes(x = Pregnant.Or.Atresia, y = pH)) + 
  geom_boxplot() +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "pH", x = "Parturition Group", y = "pH", color = "Category") +
  scale_y_continuous() + 
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))
print(pH_primary_lm_boxplot1)

## `geom_smooth()` using formula = 'y ~ x'

pH_primary_lm_boxplot2 <- ggplot(data = Primary_Samples, aes(x = Ambient.Or.OAH, y = pH)) + 
  geom_boxplot() +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "pH", x = "Treatment Group", y = "pH", color = "Category") +
  scale_y_continuous() + 
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))
print(pH_primary_lm_boxplot2)

## `geom_smooth()` using formula = 'y ~ x'

pH_primary_lm_boxplot_combo <- ggplot(data = Primary_Samples, aes(x = Pregnant.Or.Atresia, y = pH, fill = Ambient.Or.OAH)) +
  geom_boxplot(aes(x = Pregnant.Or.Atresia, y = pH, fill = Ambient.Or.OAH)) +
  geom_point(aes(color = Ambient.Or.OAH), position = position_dodge(width = 0.75), alpha = 0.3, color = "black") +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Primary pH", x = "Parturition Type", y = "pH") +
  guides(fill = guide_legend((title = "Treatment Type"))) + 
  scale_fill_manual(values = c("Ambient" = "#189bff", "OAH" = "#598c78")) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))
  
print(pH_primary_lm_boxplot_combo)

## `geom_smooth()` using formula = 'y ~ x'

## Warning: The following aesthetics were dropped during statistical transformation: fill.
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
##   the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?

# pH
# lm Boxplot

# Include both regression lines between parturition & treatment groups

ggplot(data = Primary_Samples) +
  geom_boxplot(aes(x = Pregnant.Or.Atresia, y = pH, fill = Pregnant.Or.Atresia)) +
  geom_point(aes(x = Pregnant.Or.Atresia, y = pH), position = position_dodge(width = 0.75), alpha = 0.3, color = "black") +
  geom_boxplot(aes(x = Ambient.Or.OAH, y = pH, fill = Ambient.Or.OAH)) +
  geom_point(aes(x = Ambient.Or.OAH, y = pH), position = position_dodge(width = 0.75), alpha = 0.3, color = "black")

ggplot(data = Primary_Samples) +
  geom_boxplot(aes(x = Pregnant.Or.Atresia, y = pH, fill = Pregnant.Or.Atresia)) +
  geom_point(aes(x = Pregnant.Or.Atresia, y = pH), position = position_dodge(width = 0.75), alpha = 0.3, color = "black") +
  geom_boxplot(aes(x = Ambient.Or.OAH, y = pH, fill = Ambient.Or.OAH)) +
  geom_point(aes(x = Ambient.Or.OAH, y = pH), position = position_dodge(width = 0.75), alpha = 0.3, color = "black") +
  geom_smooth(aes(x = Pregnant.Or.Atresia, y = pH, fill = Pregnant.Or.Atresia, group = 1), method = "lm") +
  geom_smooth(aes(x = Ambient.Or.OAH, y = pH, fill = Ambient.Or.OAH, group = 1), method = "lm")

## `geom_smooth()` using formula = 'y ~ x'

## Warning: The following aesthetics were dropped during statistical transformation: fill.
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
##   the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?

## `geom_smooth()` using formula = 'y ~ x'

## Warning: The following aesthetics were dropped during statistical transformation: fill.
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
##   the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?

ggplot(data = Primary_Samples) +
  geom_boxplot(aes(x = Pregnant.Or.Atresia, y = pH, fill = Pregnant.Or.Atresia)) +
  geom_point(aes(x = Pregnant.Or.Atresia, y = pH), position = position_dodge(width = 0.75), alpha = 0.3, color = "black") +
  geom_boxplot(aes(x = Ambient.Or.OAH, y = pH, fill = Ambient.Or.OAH)) +
  geom_point(aes(x = Ambient.Or.OAH, y = pH), position = position_dodge(width = 0.75), alpha = 0.3, color = "black") +
  geom_smooth(aes(x = Pregnant.Or.Atresia, y = pH, fill = Pregnant.Or.Atresia, group = 1), method = "lm") +
  geom_smooth(aes(x = Ambient.Or.OAH, y = pH, fill = Ambient.Or.OAH, group = 1), method = "lm") +
  labs(title = "Primary pH", x = "Parturition Type", y = "pH") +
  guides(fill = guide_legend((title = "Treatment Type"))) +
  scale_fill_manual(values = c("Ambient" = "#189bff", "OAH" = "#598c78")) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

## `geom_smooth()` using formula = 'y ~ x'

## Warning: The following aesthetics were dropped during statistical transformation: fill.
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
##   the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?

## `geom_smooth()` using formula = 'y ~ x'

## Warning: The following aesthetics were dropped during statistical transformation: fill.
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
##   the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?

ggplot(data = Primary_Samples) +
  geom_boxplot(aes(x = Pregnant.Or.Atresia, y = pH, fill = Ambient.Or.OAH)) +
  geom_point(aes(x = Pregnant.Or.Atresia, y = pH, fill = Ambient.Or.OAH), position = position_dodge(width = 0.75), alpha = 0.3, color = "black") +
  geom_boxplot(aes(x = Ambient.Or.OAH, y = pH, fill = Pregnant.Or.Atresia)) +
  geom_point(aes(x = Ambient.Or.OAH, y = pH, fill = Pregnant.Or.Atresia), position = position_dodge(width = 0.75), alpha = 0.3, color = "black") +
  geom_smooth(aes(x = Pregnant.Or.Atresia, y = pH, group = 1), method = "lm") +
  geom_smooth(aes(x = Ambient.Or.OAH, y = pH, group = 1), method = "lm") +
  labs(title = "Primary pH", x = "Parturition Type", y = "pH") +
  guides(fill = guide_legend((title = "Treatment Type"))) +
  scale_fill_manual(values = c("Ambient" = "#189bff", "OAH" = "#598c78")) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

Hct : Hematocrit

Hct Summary Stats

# hct
# Summary stats

# Primary data
Primary_Samples %>%
  group_by(Pregnant.Or.Atresia, Ambient.Or.OAH) %>%
  summarize(count = n(),
            median = round(median(Hct....), 3),
            mean = round(mean(Hct....), 3),
            sd = round(sd(Hct....), 3),
            cv = round(sd(Hct....)/mean(Hct....), 3)) %>%
  ungroup()

## `summarise()` has grouped output by 'Pregnant.Or.Atresia'. You can override
## using the `.groups` argument.

## # A tibble: 4 × 7
##   Pregnant.Or.Atresia Ambient.Or.OAH count median  mean    sd    cv
##   <ord>               <ord>          <int>  <dbl> <dbl> <dbl> <dbl>
## 1 Post Parturition    Ambient            8   14.3  15    3.70 0.247
## 2 Post Parturition    OAH                2   16.3  16.3  2.69 0.165
## 3 Atresia             Ambient           13   17.3  18.2  2.33 0.128
## 4 Atresia             OAH                5   18.6  20    4.50 0.225

Hct Boxplot

# hct boxplot

# Primary samples
hct_primary_boxplot <- ggplot(data = Primary_Samples, aes(x = Pregnant.Or.Atresia, y = Hct...., fill = Ambient.Or.OAH)) +
  geom_boxplot(aes(x = Pregnant.Or.Atresia, y = Hct...., fill = Ambient.Or.OAH)) +
  geom_point(aes(color = Ambient.Or.OAH), position = position_dodge(width = 0.75), alpha = 0.3, color = "black") +
  labs(title = "Hematocrit", x = "Parturition Type", y = "Hematocrit (%)") +
  guides(fill = guide_legend((title = "Treatment Type"))) + 
  scale_fill_manual(values = c("Ambient" = "#189bff", "OAH" = "#598c78")) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(hct_primary_boxplot)

ggsave(filename = "data-images/hct_primary_boxplot.pdf", plot = hct_primary_boxplot, device = "pdf")

## Saving 7 x 5 in image

ggsave(filename = "data-images/hct_primary_boxplot.png", plot = hct_primary_boxplot, device = "png")

## Saving 7 x 5 in image

Data Distribution

# hct
# Parametric Assumptions

# Data Distribution
hist(Primary_Samples$Hct....)

plotNormalHistogram(Primary_Samples$Hct....) 

plotNormalDensity(Primary_Samples$Hct....)

Hct Stat Tests

Differences:

# hct
# Primary Samples

# Interactive aov model
hct_primary_aov_int <- aov(Hct.... ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)
summary(hct_primary_aov_int) # interaction almost significant

##                                    Df Sum Sq Mean Sq F value Pr(>F)  
## Pregnant.Or.Atresia                 1  76.07   76.07   7.334 0.0123 *
## Ambient.Or.OAH                      1  14.13   14.13   1.362 0.2547  
## Pregnant.Or.Atresia:Ambient.Or.OAH  1   0.28    0.28   0.027 0.8715  
## Residuals                          24 248.94   10.37                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Additive aov model
# If interaction not significant (p-value > 0.05), use additive model
hct_primary_aov_additive <- aov(Hct.... ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = Primary_Samples)
summary(hct_primary_aov_additive)

##                     Df Sum Sq Mean Sq F value Pr(>F)  
## Pregnant.Or.Atresia  1  76.07   76.07   7.631 0.0106 *
## Ambient.Or.OAH       1  14.13   14.13   1.417 0.2451  
## Residuals           25 249.22    9.97                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Parametric Post-Hoc: Tukey's HSD test
TukeyHSD(hct_primary_aov_additive)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Hct.... ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = Primary_Samples)
## 
## $Pregnant.Or.Atresia
##                          diff       lwr      upr     p adj
## Atresia-Post Parturition 3.44 0.8753289 6.004671 0.0106014
## 
## $Ambient.Or.OAH
##                 diff       lwr      upr     p adj
## OAH-Ambient 1.634286 -1.203694 4.472265 0.2467681

# Normality plot

# Qqplot: Residuals of aov or lm model (plot residuals of aov vs lm are identical)
hct_primary_res_qqplot <- ggqqplot(residuals(aov(Hct.... ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples))) +
labs(title = "Primary Hematocrit Interactive Residual QQplot",
     subtitle = "residuals(aov(Hct.... ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples))",
                               x = "Hct Theoretical Quantiles (Predicted values)",
                               y = "Hct Sample Quantiles") +
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5))
print(hct_primary_res_qqplot)

# Normality Test
shapiro.test(Primary_Samples$Hct....)

## 
##  Shapiro-Wilk normality test
## 
## data:  Primary_Samples$Hct....
## W = 0.9866, p-value = 0.9687

shapiro.test(residuals(lm(Hct.... ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(lm(Hct.... ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples))
## W = 0.96395, p-value = 0.4305

# Parametric variance test:
bartlett.test(Hct.... ~ Pregnant.Or.Atresia, data = Primary_Samples) 

## 
##  Bartlett test of homogeneity of variances
## 
## data:  Hct.... by Pregnant.Or.Atresia
## Bartlett's K-squared = 0.16489, df = 1, p-value = 0.6847

bartlett.test(Hct.... ~ Ambient.Or.OAH, data = Primary_Samples)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  Hct.... by Ambient.Or.OAH
## Bartlett's K-squared = 0.66339, df = 1, p-value = 0.4154

# Nonparametric variance test:
# LeveneTest(Hct.... ~ Pregnant.Or.Atresia, data = Primary_Samples) 
# leveneTest(Hct.... ~ Ambient.Or.OAH, data = Primary_Samples)

# Nonparametric Stat test & Post-Hoc:

# kruskal.test(Hct.... ~ Pregnant.Or.Atresia, data = Primary_Samples) 
# kruskal.test(Hct.... ~ Ambient.Or.OAH, data = Primary_Samples) 
# 
# DunnTest(Hct.... ~ Pregnant.Or.Atresia, data = Primary_Samples)
# DunnTest(Hct.... ~ Ambient.Or.OAH, data = Primary_Samples)

Parametric Results: Ambient Fecundity No Atresia Sample

• Interaction not significant (Pr(>F) = 0.8715), additive model used.

• Additive summary results:

• A significant difference detected in pH between pregnancy groups (Pr(>F) = 0.0106 *), but not treatment groups (Pr(>F) = 0.2451)

• Parametric post Hoc: Tukey HSD with additive model

• Meaning: Atresia group displayed a significant increase in hematocrit compared to post parturition group.

• Evidence: Tukey HSD, Atresia-Post Parturition, p adj = 0.0106014

Parametric Assumptions: Across Sample comparisons…

• Sample size (n > 30): Not Met

• Normality: Met

• Homoscedasticity: Met

Non-Parametric Results: NA

Similarities:

# hct
# Primary Samples

# Regression Model 

# If interaction not significant, use additive model
# p-value significant if (a < 0.1 or a < 0.05)
hct_primary_lm_int <- lm(Hct.... ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)
summary(hct_primary_lm_int)

## 
## Call:
## lm(formula = Hct.... ~ Pregnant.Or.Atresia * Ambient.Or.OAH, 
##     data = Primary_Samples)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -5.100 -1.950 -0.900  2.025  6.700 
## 
## Coefficients:
##                                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                             17.3750     0.7647  22.723   <2e-16 ***
## Pregnant.Or.Atresia.L                    2.4395     1.0814   2.256   0.0335 *  
## Ambient.Or.OAH.L                         1.0960     1.0814   1.014   0.3209    
## Pregnant.Or.Atresia.L:Ambient.Or.OAH.L   0.2500     1.5293   0.163   0.8715    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.221 on 24 degrees of freedom
## Multiple R-squared:  0.2666, Adjusted R-squared:  0.1749 
## F-statistic: 2.908 on 3 and 24 DF,  p-value: 0.05532

hct_primary_lm_add <- lm(Hct.... ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = Primary_Samples)
summary(hct_primary_lm_add)

## 
## Call:
## lm(formula = Hct.... ~ Pregnant.Or.Atresia + Ambient.Or.OAH, 
##     data = Primary_Samples)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.0307 -2.0587 -0.9367  1.8064  6.7693 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            17.4099     0.7198   24.19   <2e-16 ***
## Pregnant.Or.Atresia.L   2.3419     0.8838    2.65   0.0138 *  
## Ambient.Or.OAH.L        1.1642     0.9780    1.19   0.2451    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.157 on 25 degrees of freedom
## Multiple R-squared:  0.2657, Adjusted R-squared:  0.207 
## F-statistic: 4.524 on 2 and 25 DF,  p-value: 0.02104

# hct_primary_lm_res_plot <- simulateResiduals(fittedModel = hct_primary_lm_int)
# plot(hct_primary_lm_res_plot)

hct_primary_lm_res_plot <- simulateResiduals(fittedModel = hct_primary_lm_add)
plot(hct_primary_lm_res_plot)

visualize(hct_primary_lm_add, plot = "model")

# Residuals of model for shapiro test
# hct_primary_lm_int_res <- residuals(hct_primary_lm_int) # does not work
hct_primary_lm_add_res <- residuals(hct_primary_lm_add)

# Normality test
#shapiro.test(hct_primary_lm_int_res) # does not work
shapiro.test(hct_primary_lm_add_res)

## 
##  Shapiro-Wilk normality test
## 
## data:  hct_primary_lm_add_res
## W = 0.96235, p-value = 0.3959

# Scedasticity test: 
# bartlett's test parametric
bartlett.test(Hct.... ~ Pregnant.Or.Atresia, data = Primary_Samples) 

## 
##  Bartlett test of homogeneity of variances
## 
## data:  Hct.... by Pregnant.Or.Atresia
## Bartlett's K-squared = 0.16489, df = 1, p-value = 0.6847

bartlett.test(Hct.... ~ Ambient.Or.OAH, data = Primary_Samples)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  Hct.... by Ambient.Or.OAH
## Bartlett's K-squared = 0.66339, df = 1, p-value = 0.4154

# Levene's test if nonparametric
# leveneTest(Hct.... ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)
# leveneTest(Hct.... ~ Pregnant.Or.Atresia, data = Primary_Samples) 
# leveneTest(Hct.... ~ Ambient.Or.OAH, data = Primary_Samples)

# Non-parametric regression test???:

Results: Primary Sample (With both treatments and parturition groups)

• lm interactive model not significant (Pr(>|t|) = 0.8715), additive model used

Linear Model: Additive lm param ~ parturition + treatment - A significant relationship detected - Evidence: Multiple R-squared: 0.2657, Adjusted R-squared: 0.207 F-statistic: 4.524 on 2 and 25 DF, p-value: 0.02104

• How to know which groups are correlated, and how strongly?
• (test lm individually, next code block)

Parametric Assumptions: Across Sample comparisons…

• Sample size (n > 30): Not Met

• Normality (of residuals) : Met

• Homoscedasticity: Met

# hct
# Primary Samples

# Regression Model pt 2

# Individually assess parameter to parturition groups & treatment groups
hct_primary_lm1 <- lm(Hct.... ~ Pregnant.Or.Atresia, data = Primary_Samples)

hct_primary_lm2 <- lm(Hct.... ~ Ambient.Or.OAH, data = Primary_Samples)

# p-value significant if (a < 0.1 or a < 0.05)
summary(hct_primary_lm1) # Significant

## 
## Call:
## lm(formula = Hct.... ~ Pregnant.Or.Atresia, data = Primary_Samples)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -5.360 -2.025 -0.980  2.075  6.600 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            16.9800     0.6276  27.055   <2e-16 ***
## Pregnant.Or.Atresia.L   2.4324     0.8876   2.741   0.0109 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.183 on 26 degrees of freedom
## Multiple R-squared:  0.2241, Adjusted R-squared:  0.1943 
## F-statistic: 7.511 on 1 and 26 DF,  p-value: 0.01094

summary(hct_primary_lm2) # not significant

## 
## Call:
## lm(formula = Hct.... ~ Ambient.Or.OAH, data = Primary_Samples)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.0810 -2.4810 -0.3119  2.6690  6.3571 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       17.9619     0.7646  23.492   <2e-16 ***
## Ambient.Or.OAH.L   1.3873     1.0813   1.283    0.211    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.504 on 26 degrees of freedom
## Multiple R-squared:  0.05954,    Adjusted R-squared:  0.02336 
## F-statistic: 1.646 on 1 and 26 DF,  p-value: 0.2108

hct_primary_lm_res_plot1 <- simulateResiduals(fittedModel = hct_primary_lm1)
plot(hct_primary_lm_res_plot1)

hct_primary_lm_res_plot2 <- simulateResiduals(fittedModel = hct_primary_lm2)
plot(hct_primary_lm_res_plot2)

visualize(hct_primary_lm1, plot = "model") 

visualize(hct_primary_lm2, plot = "model")

# For shapiro test
hct_primary_lm_res1 <- residuals(hct_primary_lm1)

hct_primary_lm_res2 <- residuals(hct_primary_lm2)

# Normality test
shapiro.test(hct_primary_lm_res1) # Normal

## 
##  Shapiro-Wilk normality test
## 
## data:  hct_primary_lm_res1
## W = 0.95107, p-value = 0.2108

shapiro.test(hct_primary_lm_res2) # Normal

## 
##  Shapiro-Wilk normality test
## 
## data:  hct_primary_lm_res2
## W = 0.97352, p-value = 0.6773

Results: Primary Sample Regression Analysis pt2

Linear Model: Individual test correlation between param and treatment, and param and parturition outcome

Correlation: Hct and Post-parturition vs atresia - A significant relationship detected with hematocrit and parturition groups - Evidence: Residual standard error: 3.183 on 26 degrees of freedom Multiple R-squared: 0.2241, Adjusted R-squared: 0.1943 F-statistic: 7.511 on 1 and 26 DF, p-value: 0.01094

• Meaning: Atresia group displays a significant correlative increase in hematocrit (p-value: 0.01094).

Warning: DHARMA states no significant (n.s.) within group deviation from uniformity, and the Levene Test for homogeneity of variance n.s.

Parametric Assumptions: Across Sample comparisons…

• Sample size (n > 30): Not Met

• Normality: Met

• Homoscedasticity: Met (check previous code chunk)

Residuals Plot (if significant)

# pH
# Primary samples

# ggplot with geom_smooth(aes(x=,y=)), method = "lm")


# lm for hct & parturition scatterplot: 
hct_primary_lm_scatterplot1 <- ggplot(data = Primary_Samples, aes(x = Pregnant.Or.Atresia, y = Hct....)) + 
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Hematocrit", x = "Parturition Group", y = "Hematocrit (%)", color = "Category") +
  scale_y_continuous() + 
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(hct_primary_lm_scatterplot1)

## `geom_smooth()` using formula = 'y ~ x'

ggsave(filename = "data-images/hct_primary_lm_scatterplot1.pdf", plot = hct_primary_lm_scatterplot1, device = "pdf")

## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'

ggsave(filename = "data-images/hct_primary_lm_scatterplot1.png", plot = hct_primary_lm_scatterplot1, device = "png")

## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'

hct_primary_lm_scatterplot2 <- ggplot(data = Primary_Samples, aes(x = Ambient.Or.OAH, y = Hct....)) + 
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Hematocrit", x = "Parturition Group", y = "Hematocrit (%)", color = "Category") +
  scale_y_continuous() + 
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(hct_primary_lm_scatterplot2)

## `geom_smooth()` using formula = 'y ~ x'

# pH
# lm Boxplot

hct_primary_lm_boxplot1 <- ggplot(data = Primary_Samples, aes(x = Pregnant.Or.Atresia, y = Hct....)) + 
  geom_boxplot() +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Hematocrit", x = "Parturition Group", y = "Hematocrit (%)", color = "Category") +
  scale_y_continuous() + 
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))
print(hct_primary_lm_boxplot1)

## `geom_smooth()` using formula = 'y ~ x'

hct_primary_lm_boxplot2 <- ggplot(data = Primary_Samples, aes(x = Ambient.Or.OAH, y = Hct....)) + 
  geom_boxplot() +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Hematocrit", x = "Parturition Group", y = "Hematocrit (%)", color = "Category") +
  scale_y_continuous() + 
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))
print(hct_primary_lm_boxplot2)

## `geom_smooth()` using formula = 'y ~ x'

Glucose

Glucose Summary Stats

# glucose
# Summary stats

# Primary data
Primary_Samples %>%
  group_by(Pregnant.Or.Atresia, Ambient.Or.OAH) %>%
  summarize(count = n(),
            median = round(median(Glu..mmol.L.), 3),
            mean = round(mean(Glu..mmol.L.), 3),
            sd = round(sd(Glu..mmol.L.), 3),
            cv = round(sd(Glu..mmol.L.)/mean(Glu..mmol.L.), 3)) %>%
  ungroup()

## `summarise()` has grouped output by 'Pregnant.Or.Atresia'. You can override
## using the `.groups` argument.

## # A tibble: 4 × 7
##   Pregnant.Or.Atresia Ambient.Or.OAH count median  mean    sd    cv
##   <ord>               <ord>          <int>  <dbl> <dbl> <dbl> <dbl>
## 1 Post Parturition    Ambient            8   1.55  1.61 0.591 0.367
## 2 Post Parturition    OAH                2   2.75  2.75 1.63  0.591
## 3 Atresia             Ambient           13   2     2.1  0.698 0.332
## 4 Atresia             OAH                5   2.1   1.98 0.356 0.18

Glucose Boxplot

# glucose boxplot

# Primary samples
glucose_primary_boxplot <- ggplot(data = Primary_Samples, aes(x = Pregnant.Or.Atresia, y = Glu..mmol.L., fill = Ambient.Or.OAH)) +
  geom_boxplot(aes(x = Pregnant.Or.Atresia, y = Glu..mmol.L., fill = Ambient.Or.OAH)) +
  geom_point(aes(color = Ambient.Or.OAH), position = position_dodge(width = 0.75), alpha = 0.3, color = "black") +
  labs(title = "Glucose", x = "Parturition Type", y = "Glucose (mmol/L)") +
  guides(fill = guide_legend((title = "Treatment Type"))) + 
  scale_fill_manual(values = c("Ambient" = "#189bff", "OAH" = "#598c78")) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(glucose_primary_boxplot)

ggsave(filename = "data-images/glucose_primary_boxplot.pdf", plot = glucose_primary_boxplot, device = "pdf")

## Saving 7 x 5 in image

ggsave(filename = "data-images/glucose_primary_boxplot.png", plot = glucose_primary_boxplot, device = "png")

## Saving 7 x 5 in image

Data Distribution

# glucose
# Parametric Assumptions

# Data Distribution
hist(Primary_Samples$Glu..mmol.L.)

plotNormalHistogram(Primary_Samples$Glu..mmol.L.) 

plotNormalDensity(Primary_Samples$Glu..mmol.L.)

Glucose Stat Tests

Differences:

# glucose
# Primary Samples

# Interactive aov model
glucose_primary_aov_int <- aov(Glu..mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)
summary(glucose_primary_aov_int) # interaction significant

##                                    Df Sum Sq Mean Sq F value Pr(>F)  
## Pregnant.Or.Atresia                 1  0.330  0.3303   0.693 0.4134  
## Ambient.Or.OAH                      1  0.369  0.3690   0.774 0.3877  
## Pregnant.Or.Atresia:Ambient.Or.OAH  1  1.753  1.7533   3.678 0.0671 .
## Residuals                          24 11.442  0.4767                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Additive aov model
# If interaction not significant (p-value > 0.05), use additive model
# glucose_primary_aov_additive <- aov(Glu..mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = Primary_Samples)
# summary(glucose_primary_aov_additive)

# Parametric Post-Hoc: Tukey's HSD test
TukeyHSD(glucose_primary_aov_int)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Glu..mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)
## 
## $Pregnant.Or.Atresia
##                               diff        lwr       upr     p adj
## Atresia-Post Parturition 0.2266667 -0.3353791 0.7887124 0.4134193
## 
## $Ambient.Or.OAH
##                 diff        lwr       upr     p adj
## OAH-Ambient 0.264127 -0.3578141 0.8860681 0.3894475
## 
## $`Pregnant.Or.Atresia:Ambient.Or.OAH`
##                                                  diff        lwr       upr
## Atresia:Ambient-Post Parturition:Ambient       0.4875 -0.3684013 1.3434013
## Post Parturition:OAH-Post Parturition:Ambient  1.1375 -0.3683119 2.6433119
## Atresia:OAH-Post Parturition:Ambient           0.3675 -0.7183564 1.4533564
## Post Parturition:OAH-Atresia:Ambient           0.6500 -0.7967373 2.0967373
## Atresia:OAH-Atresia:Ambient                   -0.1200 -1.1223290 0.8823290
## Atresia:OAH-Post Parturition:OAH              -0.7700 -2.3636015 0.8236015
##                                                   p adj
## Atresia:Ambient-Post Parturition:Ambient      0.4131342
## Post Parturition:OAH-Post Parturition:Ambient 0.1869587
## Atresia:OAH-Post Parturition:Ambient          0.7872361
## Post Parturition:OAH-Atresia:Ambient          0.6087444
## Atresia:OAH-Atresia:Ambient                   0.9872576
## Atresia:OAH-Post Parturition:OAH              0.5518182

# Normality plot

# Qqplot: Residuals of aov or lm model (plot residuals of aov vs lm are identical)
glucose_primary_res_qqplot <- ggqqplot(residuals(aov(Glu..mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples))) +
labs(title = "Glucose Interactive Residual QQplot",
     subtitle = "residuals(aov(Glu..mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples))",
                               x = "Glucose Theoretical Quantiles (Predicted values)",
                               y = "Glucose Sample Quantiles") +
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5))
print(glucose_primary_res_qqplot)

# Normality Test
shapiro.test(Primary_Samples$Glu..mmol.L.)

## 
##  Shapiro-Wilk normality test
## 
## data:  Primary_Samples$Glu..mmol.L.
## W = 0.84841, p-value = 0.0008683

shapiro.test(residuals(lm(Glu..mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(lm(Glu..mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples))
## W = 0.88835, p-value = 0.00617

# Parametric variance test:
bartlett.test(Glu..mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples) 

## 
##  Bartlett test of homogeneity of variances
## 
## data:  Glu..mmol.L. by Pregnant.Or.Atresia
## Bartlett's K-squared = 1.6702, df = 1, p-value = 0.1962

bartlett.test(Glu..mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  Glu..mmol.L. by Ambient.Or.OAH
## Bartlett's K-squared = 0.27111, df = 1, p-value = 0.6026

# Nonparametric variance test:
LeveneTest(Glu..mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.9338 0.3428
##       26

leveneTest(Glu..mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.0651 0.8005
##       26

# Nonparametric Stat test & Post-Hoc:

kruskal.test(Glu..mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  Glu..mmol.L. by Pregnant.Or.Atresia
## Kruskal-Wallis chi-squared = 2.9147, df = 1, p-value = 0.08778

kruskal.test(Glu..mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  Glu..mmol.L. by Ambient.Or.OAH
## Kruskal-Wallis chi-squared = 0.91756, df = 1, p-value = 0.3381

DunnTest(Glu..mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)

## 
##  Dunn's test of multiple comparisons using rank sums : holm  
## 
##                          mean.rank.diff   pval    
## Atresia-Post Parturition       5.522222 0.0878 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

DunnTest(Glu..mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)

## 
##  Dunn's test of multiple comparisons using rank sums : holm  
## 
##             mean.rank.diff   pval    
## OAH-Ambient       3.428571 0.3381    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Parametric Results: Ambient Fecundity No Atresia Sample

• Interaction significant (Pr(>F) = 0.0671), interactive model used.

• Interactive summary results:

• ask?

• Parametric post Hoc: Tukey HSD with interactive model

• Meaning: No groups appear significantly different from eachother…

• Evidence:

Parametric Assumptions: Across Sample comparisons…

• Sample size (n > 30): Not Met

• Normality: Not met, W = 0.84841, p-value = 0.0008683

• Homoscedasticity: Met

Non-Parametric Results:

• nonparametric scedasticity test (Levene’s test): Met

Kruskal-Wallis rank sum test: - A significant difference in glucose detected between pregnant and atresia groups - Evidence: Kruskal-Wallis chi-squared = 2.9147, df = 1, p-value = 0.08778

Nonparametric Post-Hoc: Dunn’s test - A significant increase in glucose detected in Atresia group compared to post parturition group (mean.rank.diff = 5.522222, pval = 0.0878).

Similarities:

# glucose
# Primary Samples

# Regression Model 

# If interaction not significant, use additive model
# p-value significant if (a < 0.1 or a < 0.05)
glucose_primary_lm_int <- lm(Glu..mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples) 
summary(glucose_primary_lm_int) # interaction significant

## 
## Call:
## lm(formula = Glu..mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, 
##     data = Primary_Samples)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.1500 -0.3344 -0.0400  0.2050  2.1000 
## 
## Coefficients:
##                                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                             2.11062    0.16393  12.875 2.87e-12 ***
## Pregnant.Or.Atresia.L                  -0.09988    0.23184  -0.431   0.6704    
## Ambient.Or.OAH.L                        0.35974    0.23184   1.552   0.1338    
## Pregnant.Or.Atresia.L:Ambient.Or.OAH.L -0.62875    0.32787  -1.918   0.0671 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6905 on 24 degrees of freedom
## Multiple R-squared:  0.1765, Adjusted R-squared:  0.07358 
## F-statistic: 1.715 on 3 and 24 DF,  p-value: 0.1907

glucose_primary_lm_add <- lm(Glu..mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = Primary_Samples)
summary(glucose_primary_lm_add)

## 
## Call:
## lm(formula = Glu..mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, 
##     data = Primary_Samples)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.8868 -0.4033 -0.1228  0.1087  2.2073 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             2.0228     0.1656  12.213 4.92e-12 ***
## Pregnant.Or.Atresia.L   0.1456     0.2034   0.716    0.481    
## Ambient.Or.OAH.L        0.1882     0.2250   0.836    0.411    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7265 on 25 degrees of freedom
## Multiple R-squared:  0.05033,    Adjusted R-squared:  -0.02565 
## F-statistic: 0.6624 on 2 and 25 DF,  p-value: 0.5244

# Simulate residuals plot

# glucose_primary_lm_res_plot <- simulateResiduals(fittedModel = glucose_primary_lm_int)
# plot(glucose_primary_lm_res_plot)

glucose_primary_lm_res_plot <- simulateResiduals(fittedModel = glucose_primary_lm_int)
plot(glucose_primary_lm_res_plot)

visualize(glucose_primary_lm_int, plot = "model")

# Residuals of model for shapiro test
glucose_primary_lm_int_res <- residuals(glucose_primary_lm_int) 
glucose_primary_lm_add_res <- residuals(glucose_primary_lm_add)

# Normality test
shapiro.test(glucose_primary_lm_int_res) # works

## 
##  Shapiro-Wilk normality test
## 
## data:  glucose_primary_lm_int_res
## W = 0.88835, p-value = 0.00617

shapiro.test(glucose_primary_lm_add_res) # works

## 
##  Shapiro-Wilk normality test
## 
## data:  glucose_primary_lm_add_res
## W = 0.81579, p-value = 0.0002053

# Scedasticity test: 
# bartlett's test parametric
bartlett.test(Glu..mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples) 

## 
##  Bartlett test of homogeneity of variances
## 
## data:  Glu..mmol.L. by Pregnant.Or.Atresia
## Bartlett's K-squared = 1.6702, df = 1, p-value = 0.1962

bartlett.test(Glu..mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  Glu..mmol.L. by Ambient.Or.OAH
## Bartlett's K-squared = 0.27111, df = 1, p-value = 0.6026

# Levene's test if nonparametric
leveneTest(Glu..mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3  1.9355 0.1508
##       24

leveneTest(Glu..mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.9338 0.3428
##       26

leveneTest(Glu..mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.0651 0.8005
##       26

# Non-parametric regression test???:

Results: Primary Sample (With both treatments and parturition groups)

• lm interactive model significant (Pr(>|t|) = 0.0671 .), interactive model used

Linear Model: Interactive lm param ~ parturition * treatment - No significant relationship detected across groups - Evidence: Multiple R-squared: 0.1765, Adjusted R-squared: 0.07358 F-statistic: 1.715 on 3 and 24 DF, p-value: 0.1907

• How to know which groups are correlated, and how strongly?
• (test lm individually, next code block)

Parametric Assumptions: Across Sample comparisons…

• Sample size (n > 30): Not Met

• Normality (of residuals) : Not met

• Homoscedasticity: Met

Non-parametric Tests:

• nonparametric normality test (Levene’s test): Met

• Nonparametric regression stat test???

# glucose
# Primary Samples

# Regression Model pt 2

# Individually assess parameter to parturition groups & treatment groups
glucose_primary_lm1 <- lm(Glu..mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)

glucose_primary_lm2 <- lm(Glu..mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)

# p-value significant if (a < 0.1 or a < 0.05)
summary(glucose_primary_lm1) # Significant

## 
## Call:
## lm(formula = Glu..mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.9400 -0.3850 -0.1033  0.1333  2.1333 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             1.9533     0.1424  13.714 2.06e-13 ***
## Pregnant.Or.Atresia.L   0.1603     0.2014   0.796    0.433    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7223 on 26 degrees of freedom
## Multiple R-squared:  0.02377,    Adjusted R-squared:  -0.01378 
## F-statistic: 0.6331 on 1 and 26 DF,  p-value: 0.4334

summary(glucose_primary_lm2) # not significant

## 
## Call:
## lm(formula = Glu..mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.0143 -0.4393 -0.1143  0.1857  2.2857 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        2.0571     0.1570   13.10 5.87e-13 ***
## Ambient.Or.OAH.L   0.2020     0.2221    0.91    0.371    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7197 on 26 degrees of freedom
## Multiple R-squared:  0.03085,    Adjusted R-squared:  -0.00643 
## F-statistic: 0.8275 on 1 and 26 DF,  p-value: 0.3714

glucose_primary_lm_res_plot1 <- simulateResiduals(fittedModel = glucose_primary_lm1)
plot(glucose_primary_lm_res_plot1) # Warning appears

glucose_primary_lm_res_plot2 <- simulateResiduals(fittedModel = glucose_primary_lm2)
plot(glucose_primary_lm_res_plot2)

visualize(glucose_primary_lm1, plot = "model") 

visualize(glucose_primary_lm2, plot = "model")

# For shapiro test
glucose_primary_lm_res1 <- residuals(glucose_primary_lm1)

glucose_primary_lm_res2 <- residuals(glucose_primary_lm2)

# Normality test
shapiro.test(glucose_primary_lm_res1) # Non-parametric

## 
##  Shapiro-Wilk normality test
## 
## data:  glucose_primary_lm_res1
## W = 0.79757, p-value = 9.675e-05

shapiro.test(glucose_primary_lm_res2) # Non-parametric

## 
##  Shapiro-Wilk normality test
## 
## data:  glucose_primary_lm_res2
## W = 0.8557, p-value = 0.001221

Results: Primary Sample Regression Analysis pt2

Linear Model: Individual test correlation between param and treatment, and param and parturition outcome

Correlation: - No significant relationship detected with glucose and parturition groups or treatment groups. - Evidence: - Meaning:

Warning: DHARMA states for Glu..mmol.L. ~ Pregnant.Or.Atresia QQ plot residuals, KS Test: P = 0.0271 Deviation Significant

Parametric Assumptions: Across Sample comparisons…

• Sample size (n > 30): Not Met

• Normality: Non-parametric See previous code block

• Homoscedasticity: Met see previous code block

Residuals Plot (if significant)

# glucose
# Primary samples

# ggplot with geom_smooth(aes(x=,y=)), method = "lm")


# lm for glucose & parturition scatterplot:
glucose_primary_lm_scatterplot1 <- ggplot(data = Primary_Samples, aes(x = Pregnant.Or.Atresia, y = Glu..mmol.L.)) +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Glucose", x = "Parturition Group", y = "Glucose (mmol/L)", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(glucose_primary_lm_scatterplot1)

## `geom_smooth()` using formula = 'y ~ x'

ggsave(filename = "data-images/glucose_primary_lm_scatterplot1.pdf", plot = glucose_primary_lm_scatterplot1, device = "pdf")

## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'

ggsave(filename = "data-images/glucose_primary_lm_scatterplot1.png", plot = glucose_primary_lm_scatterplot1, device = "png")

## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'

# lm for glucose & treatment scatterplot:
glucose_primary_lm_scatterplot2 <- ggplot(data = Primary_Samples, aes(x = Ambient.Or.OAH, y = Glu..mmol.L.)) +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Glucose", x = "Treatment Group", y = "Glucose (mmol/L)", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(glucose_primary_lm_scatterplot2)

## `geom_smooth()` using formula = 'y ~ x'

ggsave(filename = "data-images/glucose_primary_lm_scatterplot2.pdf", plot = glucose_primary_lm_scatterplot2, device = "pdf")

## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'

ggsave(filename = "data-images/glucose_primary_lm_scatterplot2.png", plot = glucose_primary_lm_scatterplot2, device = "png")

## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'

# glucose
# lm Boxplot

glucose_primary_lm_boxplot1 <- ggplot(data = Primary_Samples, aes(x = Pregnant.Or.Atresia, y = Glu..mmol.L.)) +
  geom_boxplot() +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Glucose", x = "Parturition Group", y = "Glucose (mmol/L)", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))
print(glucose_primary_lm_boxplot1)

## `geom_smooth()` using formula = 'y ~ x'

glucose_primary_lm_boxplot2 <- ggplot(data = Primary_Samples, aes(x = Ambient.Or.OAH, y = Glu..mmol.L.)) +
  geom_boxplot() +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Glucose", x = "Parturition Group", y = "Glucose (mmol/L)", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))
print(glucose_primary_lm_boxplot2)

## `geom_smooth()` using formula = 'y ~ x'

Na+ : Sodium

Sodium Summary Stats

# Na+
# Summary stats

# Primary data
Primary_Samples %>%
  group_by(Pregnant.Or.Atresia, Ambient.Or.OAH) %>%
  summarize(count = n(),
            median = round(median(Na...mmol.L.), 3),
            mean = round(mean(Na...mmol.L.), 3),
            sd = round(sd(Na...mmol.L.), 3),
            cv = round(sd(Na...mmol.L.)/mean(Na...mmol.L.), 3)) %>%
  ungroup()

## `summarise()` has grouped output by 'Pregnant.Or.Atresia'. You can override
## using the `.groups` argument.

## # A tibble: 4 × 7
##   Pregnant.Or.Atresia Ambient.Or.OAH count median  mean    sd    cv
##   <ord>               <ord>          <int>  <dbl> <dbl> <dbl> <dbl>
## 1 Post Parturition    Ambient            8   176.  176.  3.41 0.019
## 2 Post Parturition    OAH                2   178   178   0    0    
## 3 Atresia             Ambient           13   178   181. 11.4  0.063
## 4 Atresia             OAH                5   174   174   2.24 0.013

Sodium Boxplot

# Na+ boxplot

# Primary samples
sodium_primary_boxplot <- ggplot(data = Primary_Samples, aes(x = Pregnant.Or.Atresia, y = Na...mmol.L., fill = Ambient.Or.OAH)) +
  geom_boxplot(aes(x = Pregnant.Or.Atresia, y = Na...mmol.L., fill = Ambient.Or.OAH)) +
  geom_point(aes(color = Ambient.Or.OAH), position = position_dodge(width = 0.75), alpha = 0.3, color = "black") +
  labs(title = "Sodium", x = "Parturition Type", y = "Sodium (mmol/L)") +
  guides(fill = guide_legend((title = "Treatment Type"))) + 
  scale_fill_manual(values = c("Ambient" = "#189bff", "OAH" = "#598c78")) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(sodium_primary_boxplot)

ggsave(filename = "data-images/sodium_primary_boxplot.pdf", plot = sodium_primary_boxplot, device = "pdf")

## Saving 7 x 5 in image

ggsave(filename = "data-images/sodium_primary_boxplot.png", plot = sodium_primary_boxplot, device = "png")

## Saving 7 x 5 in image

Data Distribution

# Na+
# Parametric Assumptions

# Data Distribution
hist(Primary_Samples$Na...mmol.L.)

plotNormalHistogram(Primary_Samples$Na...mmol.L.) 

plotNormalDensity(Primary_Samples$Na...mmol.L.)

Na+ Stat Tests

Differences:

# Na+
# Primary Samples

# Interactive aov model
sodium_primary_aov_int <- aov(Na...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)
summary(sodium_primary_aov_int) # interaction not significant

##                                    Df Sum Sq Mean Sq F value Pr(>F)
## Pregnant.Or.Atresia                 1   58.7   58.72   0.849  0.366
## Ambient.Or.OAH                      1   97.2   97.24   1.406  0.247
## Pregnant.Or.Atresia:Ambient.Or.OAH  1   99.7   99.66   1.441  0.242
## Residuals                          24 1659.8   69.16

# Additive aov model
# If interaction not significant (p-value > 0.05), use additive model
sodium_primary_aov_add <- aov(Na...mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = Primary_Samples)
summary(sodium_primary_aov_add)

##                     Df Sum Sq Mean Sq F value Pr(>F)
## Pregnant.Or.Atresia  1   58.7   58.72   0.834  0.370
## Ambient.Or.OAH       1   97.2   97.24   1.382  0.251
## Residuals           25 1759.5   70.38

# Parametric Post-Hoc: Tukey's HSD test
TukeyHSD(sodium_primary_aov_add)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Na...mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = Primary_Samples)
## 
## $Pregnant.Or.Atresia
##                              diff       lwr     upr     p adj
## Atresia-Post Parturition 3.022222 -3.792266 9.83671 0.3697553
## 
## $Ambient.Or.OAH
##                  diff       lwr      upr     p adj
## OAH-Ambient -4.287831 -11.82852 3.252855 0.2525973

# Normality plot

# Qqplot: Residuals of aov or lm model (plot residuals of aov vs lm are identical)
sodium_primary_res_qqplot <- ggqqplot(residuals(aov(Na...mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = Primary_Samples))) +
labs(title = "Sodium Interactive Residual QQplot",
     subtitle = "residuals(aov(Na...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples))",
                               x = "Sodium Theoretical Quantiles (Predicted values)",
                               y = "Sodium Sample Quantiles") +
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5))
print(sodium_primary_res_qqplot)

# Normality Test
shapiro.test(Primary_Samples$Na...mmol.L.)

## 
##  Shapiro-Wilk normality test
## 
## data:  Primary_Samples$Na...mmol.L.
## W = 0.52152, p-value = 1.916e-08

shapiro.test(residuals(lm(Na...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(lm(Na...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples))
## W = 0.56167, p-value = 5.246e-08

# Parametric variance test:
bartlett.test(Na...mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples) 

## 
##  Bartlett test of homogeneity of variances
## 
## data:  Na...mmol.L. by Pregnant.Or.Atresia
## Bartlett's K-squared = 10.874, df = 1, p-value = 0.0009751

bartlett.test(Na...mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  Na...mmol.L. by Ambient.Or.OAH
## Bartlett's K-squared = 8.4559, df = 1, p-value = 0.003639

# Nonparametric variance test:
LeveneTest(Na...mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.4826 0.4934
##       26

leveneTest(Na...mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.2851 0.5979
##       26

# Nonparametric Stat test & Post-Hoc:

kruskal.test(Na...mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  Na...mmol.L. by Pregnant.Or.Atresia
## Kruskal-Wallis chi-squared = 0.084351, df = 1, p-value = 0.7715

kruskal.test(Na...mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  Na...mmol.L. by Ambient.Or.OAH
## Kruskal-Wallis chi-squared = 2.8469, df = 1, p-value = 0.09155

DunnTest(Na...mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)

## 
##  Dunn's test of multiple comparisons using rank sums : holm  
## 
##                          mean.rank.diff   pval    
## Atresia-Post Parturition      0.9333333 0.7715    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

DunnTest(Na...mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)

## 
##  Dunn's test of multiple comparisons using rank sums : holm  
## 
##             mean.rank.diff   pval    
## OAH-Ambient             -6 0.0916 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

sodium_primary_aov_res_plot <- simulateResiduals(fittedModel = sodium_primary_aov_add)

## Warning in checkModel(fittedModel): DHARMa: fittedModel not in class of
## supported models. Absolutely no guarantee that this will work!

plot(sodium_primary_aov_res_plot) # KS significant

visualize(sodium_primary_aov_add)

boxplot(data = Primary_Samples, Na...mmol.L. ~ Ambient.Or.OAH)

# Test outliers
outlierTest(sodium_primary_aov_add) 

##    rstudent unadjusted p-value Bonferroni p
## 24 12.31196         7.3338e-12   2.0535e-10

#  Used Ranked model to adjust for outlier
boxplot(rank(Na...mmol.L.) ~ Pregnant.Or.Atresia, data = Primary_Samples, horizontal = TRUE)

boxplot(rank(Na...mmol.L.) ~ Ambient.Or.OAH, data = Primary_Samples, horizontal = TRUE)

# Nonparametric Wilcoxon test
wilcox.test(rank(Na...mmol.L.) ~ Pregnant.Or.Atresia, data = Primary_Samples)

## Warning in wilcox.test.default(x = DATA[[1L]], y = DATA[[2L]], ...): cannot
## compute exact p-value with ties

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  rank(Na...mmol.L.) by Pregnant.Or.Atresia
## W = 84, p-value = 0.7901
## alternative hypothesis: true location shift is not equal to 0

wilcox.test(rank(Na...mmol.L.) ~ Ambient.Or.OAH, data = Primary_Samples)

## Warning in wilcox.test.default(x = DATA[[1L]], y = DATA[[2L]], ...): cannot
## compute exact p-value with ties

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  rank(Na...mmol.L.) by Ambient.Or.OAH
## W = 105, p-value = 0.09682
## alternative hypothesis: true location shift is not equal to 0

Parametric Results: Ambient Fecundity No Atresia Sample

Warning: Severe outlier may be interfering with results

• Interaction not significant (Pr(>F) = 0.242), additive model used.

• Additive summary results:

• No significant differences between parturition and treatment groups.

• Parametric post Hoc: Tukey HSD with interactive model

• Meaning: No groups appear significantly different from eachother…

• Evidence:

Parametric Assumptions: Across Sample comparisons…

• Sample size (n > 30): Not Met

• Normality: Really not met, W = 0.52152, p-value = 1.916e-08

• Homoscedasticity: Not met (for both factor groups)

Non-Parametric Results:

• nonparametric scedasticity test (Levene’s test): Met (for both groups)

Kruskal-Wallis rank sum test: - A significant difference in sodium detected between Treatments: Ambient and OAH groups - Evidence: Kruskal-Wallis chi-squared = 2.8469, df = 1, p-value = 0.09155

Nonparametric Post-Hoc: Dunn’s test - A significant difference in sodium detected in treatment groups (mean.rank.diff = -6 pval = 0.0916). - How to know which group is greater? - Generate boxplot for sodium just with treatment groups?

Similarities:

# sodium
# Primary Samples

# Regression Model 
# p-value significant if (a < 0.1 or a < 0.05)

sodium_primary_lm_int <- lm(Na...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples) # not significant
summary(sodium_primary_lm_int)

## 
## Call:
## lm(formula = Na...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, 
##     data = Primary_Samples)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -7.231 -3.481 -1.115  1.063 36.769 
## 
## Coefficients:
##                                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                            177.2452     1.9745  89.769   <2e-16 ***
## Pregnant.Or.Atresia.L                    0.5235     2.7923   0.187    0.853    
## Ambient.Or.OAH.L                        -1.7610     2.7923  -0.631    0.534    
## Pregnant.Or.Atresia.L:Ambient.Or.OAH.L  -4.7404     3.9489  -1.200    0.242    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.316 on 24 degrees of freedom
## Multiple R-squared:  0.1335, Adjusted R-squared:  0.02514 
## F-statistic: 1.232 on 3 and 24 DF,  p-value: 0.3198

sodium_primary_lm_add <- lm(Na...mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = Primary_Samples)
summary(sodium_primary_lm_add)

## 
## Call:
## lm(formula = Na...mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, 
##     data = Primary_Samples)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -8.064 -3.422 -1.743  0.907 37.578 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            176.583      1.913  92.330   <2e-16 ***
## Pregnant.Or.Atresia.L    2.375      2.348   1.011    0.322    
## Ambient.Or.OAH.L        -3.055      2.599  -1.175    0.251    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.389 on 25 degrees of freedom
## Multiple R-squared:  0.08142,    Adjusted R-squared:  0.007938 
## F-statistic: 1.108 on 2 and 25 DF,  p-value: 0.3459

# sodium_primary_lm_res_plot <- simulateResiduals(fittedModel = sodium_primary_lm_int)
# plot(sodium_primary_lm_res_plot)

sodium_primary_lm_res_plot <- simulateResiduals(fittedModel = sodium_primary_lm_add)
plot(sodium_primary_lm_res_plot) # KS significant

visualize(sodium_primary_lm_add, plot = "model")

# Check Dispersion or Outliers
plotResiduals(sodium_primary_lm_add)

testDispersion(sodium_primary_lm_add)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.9361, p-value = 0.832
## alternative hypothesis: two.sided

# Transform data: square root
sodium_primary_sqrt <- Primary_Samples %>%
  mutate(Na...mmol.L. = sqrt(Na...mmol.L.))

# Rerun lm model with square root transformed data
sodium_primary_sqrt_lm_int <- lm(Na...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = sodium_primary_sqrt) 
summary(sodium_primary_sqrt_lm_int) # still not significant

## 
## Call:
## lm(formula = Na...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, 
##     data = sodium_primary_sqrt)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.26562 -0.12425 -0.03893  0.04045  1.30830 
## 
## Coefficients:
##                                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                            13.31135    0.07073 188.202   <2e-16 ***
## Pregnant.Or.Atresia.L                   0.01733    0.10003   0.173    0.864    
## Ambient.Or.OAH.L                       -0.06389    0.10003  -0.639    0.529    
## Pregnant.Or.Atresia.L:Ambient.Or.OAH.L -0.17549    0.14146  -1.241    0.227    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2979 on 24 degrees of freedom
## Multiple R-squared:  0.1388, Adjusted R-squared:  0.03118 
## F-statistic:  1.29 on 3 and 24 DF,  p-value: 0.3006

sodium_primary_sqrt_lm_add <- lm(Na...mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = sodium_primary_sqrt)
summary(sodium_primary_sqrt_lm_add)

## 
## Call:
## lm(formula = Na...mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, 
##     data = sodium_primary_sqrt)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.30517 -0.12246 -0.06196  0.03584  1.33823 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           13.28684    0.06864 193.567   <2e-16 ***
## Pregnant.Or.Atresia.L  0.08586    0.08428   1.019    0.318    
## Ambient.Or.OAH.L      -0.11178    0.09327  -1.198    0.242    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3011 on 25 degrees of freedom
## Multiple R-squared:  0.08361,    Adjusted R-squared:  0.01029 
## F-statistic:  1.14 on 2 and 25 DF,  p-value: 0.3358

# Check hows it
sodium_primary_sqrt_lm_res_plot <- simulateResiduals(fittedModel = sodium_primary_sqrt_lm_add)
plot(sodium_primary_sqrt_lm_res_plot) # KS still significant

## qu = 0.5, log(sigma) = -3.042338 : outer Newton did not converge fully.

plotResiduals(sodium_primary_sqrt_lm_add)

## qu = 0.5, log(sigma) = -3.042338 : outer Newton did not converge fully.

testDispersion(sodium_primary_sqrt_lm_add)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.9361, p-value = 0.832
## alternative hypothesis: two.sided

# Residuals of model for shapiro test
sodium_primary_lm_int_res <- residuals(sodium_primary_lm_int) 
sodium_primary_lm_add_res <- residuals(sodium_primary_lm_add)


# Normality test
shapiro.test(sodium_primary_lm_int_res) # significant

## 
##  Shapiro-Wilk normality test
## 
## data:  sodium_primary_lm_int_res
## W = 0.56167, p-value = 5.246e-08

shapiro.test(sodium_primary_lm_add_res) # significant

## 
##  Shapiro-Wilk normality test
## 
## data:  sodium_primary_lm_add_res
## W = 0.56739, p-value = 6.085e-08

# Scedasticity test: 
# bartlett's test parametric
bartlett.test(Na...mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples) 

## 
##  Bartlett test of homogeneity of variances
## 
## data:  Na...mmol.L. by Pregnant.Or.Atresia
## Bartlett's K-squared = 10.874, df = 1, p-value = 0.0009751

bartlett.test(Na...mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  Na...mmol.L. by Ambient.Or.OAH
## Bartlett's K-squared = 8.4559, df = 1, p-value = 0.003639

# Levene's test if nonparametric
leveneTest(Na...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3  0.3841 0.7654
##       24

leveneTest(Na...mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.4826 0.4934
##       26

leveneTest(Na...mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.2851 0.5979
##       26

# Non-parametric regression test???:

Results: Primary Sample (With both treatments and parturition groups)

• lm interactive model not significant (Pr(>|t|) = 0.242), additive model used

Linear Model: Additive lm param ~ parturition + treatment (no sqrt adjust) - No significant relationship detected across groups - Evidence: Multiple R-squared: 0.08142, Adjusted R-squared: 0.007938 F-statistic: 1.108 on 2 and 25 DF, p-value: 0.3459

• How to know which groups are correlated, and how strongly?
• (test lm individually, next code block)

Warning: DHARMA states for Na…mmol.L ~ Pregnant.Or.Atresia QQ plot residuals, KS Test: P = 0.0271 Deviation Significant

Correlation: Using sqrt transformed sodium values (additive model) - No significant relationship detected with glucose and parturition groups or treatment groups. - Evidence: Multiple R-squared: 0.08361, Adjusted R-squared: 0.01029 F-statistic: 1.14 on 2 and 25 DF, p-value: 0.3358

Parametric Assumptions: Across Sample comparisons…

• Sample size (n > 30): Not Met

• Normality (of residuals) : Not met

• Homoscedasticity: Not Met

Non-parametric Tests:

• nonparametric normality test (Levene’s test): Met

• Nonparametric regression stat test???

# glucose
# Primary Samples

# Regression Model pt 2

# Individually assess parameter to parturition groups & treatment groups
sodium_primary_lm1 <- lm(Na...mmol.L. ~ Pregnant.Or.Atresia, data = sodium_primary_sqrt)

sodium_primary_lm2 <- lm(Na...mmol.L. ~ Ambient.Or.OAH, data = sodium_primary_sqrt)

# p-value significant if (a < 0.1 or a < 0.05)
summary(sodium_primary_lm1) 

## 
## Call:
## lm(formula = Na...mmol.L. ~ Pregnant.Or.Atresia, data = sodium_primary_sqrt)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.30599 -0.11729 -0.04291  0.04232  1.38214 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           13.32812    0.05987 222.606   <2e-16 ***
## Pregnant.Or.Atresia.L  0.07717    0.08467   0.911     0.37    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3036 on 26 degrees of freedom
## Multiple R-squared:  0.03096,    Adjusted R-squared:  -0.006314 
## F-statistic: 0.8306 on 1 and 26 DF,  p-value: 0.3705

summary(sodium_primary_lm2) 

## 
## Call:
## lm(formula = Na...mmol.L. ~ Ambient.Or.OAH, data = sodium_primary_sqrt)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.38034 -0.08912 -0.04079 -0.00125  1.38449 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      13.30708    0.06575 202.382   <2e-16 ***
## Ambient.Or.OAH.L -0.10360    0.09299  -1.114    0.275    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3013 on 26 degrees of freedom
## Multiple R-squared:  0.04557,    Adjusted R-squared:  0.008856 
## F-statistic: 1.241 on 1 and 26 DF,  p-value: 0.2754

sodium_primary_lm_res_plot1 <- simulateResiduals(fittedModel = sodium_primary_lm1)
plot(sodium_primary_lm_res_plot1) # Warning appears

sodium_primary_lm_res_plot2 <- simulateResiduals(fittedModel = sodium_primary_lm2)
plot(sodium_primary_lm_res_plot2)

visualize(sodium_primary_lm1, plot = "model") 

visualize(sodium_primary_lm2, plot = "model")

# For shapiro test
sodium_primary_lm_res1 <- residuals(sodium_primary_lm1)

sodium_primary_lm_res2 <- residuals(sodium_primary_lm2)

# Normality test
shapiro.test(sodium_primary_lm_res1) # Non-parametric

## 
##  Shapiro-Wilk normality test
## 
## data:  sodium_primary_lm_res1
## W = 0.58314, p-value = 9.212e-08

shapiro.test(sodium_primary_lm_res2) # Non-parametric

## 
##  Shapiro-Wilk normality test
## 
## data:  sodium_primary_lm_res2
## W = 0.55298, p-value = 4.198e-08

Results: Primary Sample Regression Analysis pt2

Linear Model: Individual test correlation between param and treatment, and param and parturition outcome

Correlation: Using sqrt transformed sodium values - No significant relationship detected with glucose and parturition groups or treatment groups. - Evidence: - Meaning:

Warning: DHARMA states for Na…mmol.L. ~ Pregnant.Or.Atresia QQ plot residuals, KS Test: P < 0.05. Deviation Significant

Parametric Assumptions: Across Sample comparisons…

• Sample size (n > 30): Not Met

• Normality: Non-parametric See previous code block

• Homoscedasticity: Not Met see previous code block

Residuals Plot (if significant)

# sodium
# Primary samples

# ggplot with geom_smooth(aes(x=,y=)), method = "lm")


# lm for sodium & parturition scatterplot:
sodium_primary_lm_scatterplot1 <- ggplot(data = sodium_primary_sqrt, aes(x = Pregnant.Or.Atresia, y = Na...mmol.L.)) +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Sodium", x = "Parturition Group", y = "Sqrt Sodium", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(sodium_primary_lm_scatterplot1)

## `geom_smooth()` using formula = 'y ~ x'

ggsave(filename = "data-images/sodium_primary_lm_scatterplot1.pdf", plot = sodium_primary_lm_scatterplot1, device = "pdf")

## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'

ggsave(filename = "data-images/sodium_primary_lm_scatterplot1.png", plot = sodium_primary_lm_scatterplot1, device = "png")

## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'

# lm for sodium & treatment scatterplot:
sodium_primary_lm_scatterplot2 <- ggplot(data = sodium_primary_sqrt, aes(x = Ambient.Or.OAH, y = Na...mmol.L.)) +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Sodium", x = "Treatment Group", y = "Sqrt Sodium", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(sodium_primary_lm_scatterplot2)

## `geom_smooth()` using formula = 'y ~ x'

ggsave(filename = "data-images/sodium_primary_lm_scatterplot2.pdf", plot = sodium_primary_lm_scatterplot2, device = "pdf")

## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'

ggsave(filename = "data-images/sodium_primary_lm_scatterplot2.png", plot = sodium_primary_lm_scatterplot2, device = "png")

## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'

# sodium
# lm Boxplot

sodium_primary_lm_boxplot1 <- ggplot(data = sodium_primary_sqrt, aes(x = Pregnant.Or.Atresia, y = Na...mmol.L.)) +
  geom_boxplot() +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Sodium", x = "Parturition Group", y = "Sqrt Sodium", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))
print(sodium_primary_lm_boxplot1)

## `geom_smooth()` using formula = 'y ~ x'

sodium_primary_lm_boxplot2 <- ggplot(data = sodium_primary_sqrt, aes(x = Ambient.Or.OAH, y = Na...mmol.L.)) +
  geom_boxplot() +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Sodium", x = "Parturition Group", y = "Sqrt Sodium", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))
print(sodium_primary_lm_boxplot2)

## `geom_smooth()` using formula = 'y ~ x'

Cl- : Chloride

Chloride Summary Stats

# Cl-
# Summary stats

# Primary data
Primary_Samples %>%
  group_by(Pregnant.Or.Atresia, Ambient.Or.OAH) %>%
  summarize(count = n(),
            median = round(median(Cl...mmol.L.), 3),
            mean = round(mean(Cl...mmol.L.), 3),
            sd = round(sd(Cl...mmol.L.), 3),
            cv = round(sd(Cl...mmol.L.)/mean(Cl...mmol.L.), 3)) %>%
  ungroup()

## `summarise()` has grouped output by 'Pregnant.Or.Atresia'. You can override
## using the `.groups` argument.

## # A tibble: 4 × 7
##   Pregnant.Or.Atresia Ambient.Or.OAH count median  mean    sd    cv
##   <ord>               <ord>          <int>  <dbl> <dbl> <dbl> <dbl>
## 1 Post Parturition    Ambient            8    157  158.  3.10 0.02 
## 2 Post Parturition    OAH                2    160  160   2.83 0.018
## 3 Atresia             Ambient           13    158  158.  4.61 0.029
## 4 Atresia             OAH                5    159  159.  1.79 0.011

Chloride Boxplot

# Cl- boxplot

# Primary samples
chloride_primary_boxplot <- ggplot(data = Primary_Samples, aes(x = Pregnant.Or.Atresia, y = Cl...mmol.L., fill = Ambient.Or.OAH)) +
  geom_boxplot(aes(x = Pregnant.Or.Atresia, y = Cl...mmol.L., fill = Ambient.Or.OAH)) +
  geom_point(aes(color = Ambient.Or.OAH), position = position_dodge(width = 0.75), alpha = 0.3, color = "black") +
  labs(title = "Chloride", x = "Parturition Type", y = "Chloride (mmol/L)") +
  guides(fill = guide_legend((title = "Treatment Type"))) + 
  scale_fill_manual(values = c("Ambient" = "#189bff", "OAH" = "#598c78")) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(chloride_primary_boxplot)

ggsave(filename = "data-images/chloride_primary_boxplot.pdf", plot = chloride_primary_boxplot, device = "pdf")

## Saving 7 x 5 in image

ggsave(filename = "data-images/chloride_primary_boxplot.png", plot = chloride_primary_boxplot, device = "png")

## Saving 7 x 5 in image

Data Distribution

# Cl-
# Parametric Assumptions

# Data Distribution
hist(Primary_Samples$Cl...mmol.L.)

plotNormalHistogram(Primary_Samples$Cl...mmol.L.) 

plotNormalDensity(Primary_Samples$Cl...mmol.L.)

Cl- Stat Tests

Differences:

# Cl-
# Primary Samples

# Interactive aov model
chloride_primary_aov_int <- aov(Cl...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)
summary(chloride_primary_aov_int) # interaction not significant

##                                    Df Sum Sq Mean Sq F value Pr(>F)
## Pregnant.Or.Atresia                 1    1.4   1.400   0.098  0.757
## Ambient.Or.OAH                      1    7.5   7.536   0.526  0.475
## Pregnant.Or.Atresia:Ambient.Or.OAH  1    2.5   2.533   0.177  0.678
## Residuals                          24  343.5  14.314

# Additive aov model
# If interaction not significant (p-value > 0.05), use additive model
chloride_primary_aov_add <- aov(Cl...mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = Primary_Samples)
summary(chloride_primary_aov_add)

##                     Df Sum Sq Mean Sq F value Pr(>F)
## Pregnant.Or.Atresia  1    1.4   1.400   0.101  0.753
## Ambient.Or.OAH       1    7.5   7.536   0.544  0.467
## Residuals           25  346.1  13.843

# Parametric Post-Hoc: Tukey's HSD test
TukeyHSD(chloride_primary_aov_add)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Cl...mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = Primary_Samples)
## 
## $Pregnant.Or.Atresia
##                               diff       lwr     upr     p adj
## Atresia-Post Parturition 0.4666667 -2.555517 3.48885 0.7531103
## 
## $Ambient.Or.OAH
##                 diff       lwr      upr     p adj
## OAH-Ambient 1.193651 -2.150597 4.537899 0.4691166

# Normality plot

# Qqplot: Residuals of aov or lm model (plot residuals of aov vs lm are identical)
chloride_primary_res_qqplot <- ggqqplot(residuals(aov(Cl...mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = Primary_Samples))) +
labs(title = "Chloride Interactive Residual QQplot",
     subtitle = "residuals(aov(Cl...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples))",
                               x = "Chloride Theoretical Quantiles (Predicted values)",
                               y = "Chloride Sample Quantiles") +
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(size = 10, hjust = 0.5))
print(chloride_primary_res_qqplot)

# Normality Test
shapiro.test(Primary_Samples$Cl...mmol.L.)

## 
##  Shapiro-Wilk normality test
## 
## data:  Primary_Samples$Cl...mmol.L.
## W = 0.97188, p-value = 0.6318

shapiro.test(residuals(lm(Cl...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(lm(Cl...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples))
## W = 0.96898, p-value = 0.5534

# Parametric variance test:
bartlett.test(Cl...mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples) 

## 
##  Bartlett test of homogeneity of variances
## 
## data:  Cl...mmol.L. by Pregnant.Or.Atresia
## Bartlett's K-squared = 0.76057, df = 1, p-value = 0.3832

bartlett.test(Cl...mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  Cl...mmol.L. by Ambient.Or.OAH
## Bartlett's K-squared = 3.6573, df = 1, p-value = 0.05582

# Nonparametric variance test:
LeveneTest(Cl...mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.1495 0.7021
##       26

leveneTest(Cl...mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value  Pr(>F)  
## group  1  2.9761 0.09637 .
##       26                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Nonparametric Stat test & Post-Hoc:

kruskal.test(Cl...mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  Cl...mmol.L. by Pregnant.Or.Atresia
## Kruskal-Wallis chi-squared = 0.23333, df = 1, p-value = 0.6291

kruskal.test(Cl...mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  Cl...mmol.L. by Ambient.Or.OAH
## Kruskal-Wallis chi-squared = 1.1429, df = 1, p-value = 0.285

DunnTest(Cl...mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)

## 
##  Dunn's test of multiple comparisons using rank sums : holm  
## 
##                          mean.rank.diff   pval    
## Atresia-Post Parturition       1.555556 0.6291    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

DunnTest(Cl...mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)

## 
##  Dunn's test of multiple comparisons using rank sums : holm  
## 
##             mean.rank.diff   pval    
## OAH-Ambient       3.809524 0.2850    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

visualize(chloride_primary_aov_add)

# Test outliers
outlierTest(chloride_primary_aov_add) 

## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##    rstudent unadjusted p-value Bonferroni p
## 23 2.316521           0.029385      0.82279

#  Used Ranked model to adjust for outlier
boxplot(rank(Cl...mmol.L.) ~ Pregnant.Or.Atresia, data = Primary_Samples, horizontal = TRUE)

boxplot(rank(Cl...mmol.L.) ~ Ambient.Or.OAH, data = Primary_Samples, horizontal = TRUE)

# Nonparametric Wilcoxon test
wilcox.test(rank(Cl...mmol.L.) ~ Pregnant.Or.Atresia, data = Primary_Samples)

## Warning in wilcox.test.default(x = DATA[[1L]], y = DATA[[2L]], ...): cannot
## compute exact p-value with ties

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  rank(Cl...mmol.L.) by Pregnant.Or.Atresia
## W = 80, p-value = 0.6463
## alternative hypothesis: true location shift is not equal to 0

wilcox.test(rank(Cl...mmol.L.) ~ Ambient.Or.OAH, data = Primary_Samples)

## Warning in wilcox.test.default(x = DATA[[1L]], y = DATA[[2L]], ...): cannot
## compute exact p-value with ties

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  rank(Cl...mmol.L.) by Ambient.Or.OAH
## W = 53.5, p-value = 0.2973
## alternative hypothesis: true location shift is not equal to 0

Parametric Results: Ambient Fecundity No Atresia Sample

• Interaction not significant (Pr(>F) = 0.678), additive model used.

• Additive summary results:

• No significant differences between parturition and treatment groups.

• Parametric post Hoc: Tukey HSD with additive model

• No groups appear significantly different from eachother…

• Evidence:

Parametric Assumptions: Across Sample comparisons…

• Sample size (n > 30): Not Met

• Normality: Met, W = 0.97188, p-value = 0.6318

• Homoscedasticity: Not met (for treatment group only)

• Evidence:

• Parturition group: Bartlett’s K-squared = 0.76057, df = 1, p-value = 0.3832

• Treatment group: Bartlett’s K-squared = 3.6573, df = 1, p-value = 0.05582

Non-Parametric Results:

• nonparametric scedasticity test (Levene’s test): Met ish
• Evidence:
• Parturition group: F value = 0.1495, Pr(>F) = 0.7021
• Treatment group: F value = 2.9761 , Pr(>F) = 0.09637

Kruskal-Wallis rank sum test: - No significant difference in chloride detected between factor groups. - Evidence:

Nonparametric Post-Hoc: Dunn’s test - No significant difference detected across factor groups. - Evidence: - Parturition group: pval = 0.6291 - Treatment group: pval = 0.2850

Similarities:

# chloride
# Primary Samples

# Regression Model 
# p-value significant if (a < 0.1 or a < 0.05)

chloride_primary_lm_int <- lm(Cl...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples) # not significant
summary(chloride_primary_lm_int)

## 
## Call:
## lm(formula = Cl...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, 
##     data = Primary_Samples)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -7.462 -2.050 -0.200  1.654  7.538 
## 
## Coefficients:
##                                         Estimate Std. Error t value Pr(>|t|)
## (Intercept)                            158.85288    0.89826 176.845   <2e-16
## Pregnant.Or.Atresia.L                   -0.03128    1.27033  -0.025    0.981
## Ambient.Or.OAH.L                         1.05658    1.27033   0.832    0.414
## Pregnant.Or.Atresia.L:Ambient.Or.OAH.L  -0.75577    1.79652  -0.421    0.678
##                                           
## (Intercept)                            ***
## Pregnant.Or.Atresia.L                     
## Ambient.Or.OAH.L                          
## Pregnant.Or.Atresia.L:Ambient.Or.OAH.L    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.783 on 24 degrees of freedom
## Multiple R-squared:  0.03231,    Adjusted R-squared:  -0.08865 
## F-statistic: 0.2671 on 3 and 24 DF,  p-value: 0.8484

chloride_primary_lm_add <- lm(Cl...mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = Primary_Samples)
summary(chloride_primary_lm_add)

## 
## Call:
## lm(formula = Cl...mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, 
##     data = Primary_Samples)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.3326 -2.0528 -0.4339  1.8667  7.6674 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           158.7473     0.8482 187.160   <2e-16 ***
## Pregnant.Or.Atresia.L   0.2638     1.0415   0.253    0.802    
## Ambient.Or.OAH.L        0.8503     1.1525   0.738    0.467    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.721 on 25 degrees of freedom
## Multiple R-squared:  0.02517,    Adjusted R-squared:  -0.05281 
## F-statistic: 0.3228 on 2 and 25 DF,  p-value: 0.7271

# chloride_primary_lm_res_plot <- simulateResiduals(fittedModel = chloride_primary_lm_int)
# plot(chloride_primary_lm_res_plot)

chloride_primary_lm_res_plot <- simulateResiduals(fittedModel = chloride_primary_lm_add)
plot(chloride_primary_lm_res_plot) # quantile deviations detected

visualize(chloride_primary_lm_add, plot = "model")

# Check Dispersion or Outliers
plotResiduals(chloride_primary_lm_add)

testDispersion(chloride_primary_lm_add)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.9361, p-value = 0.832
## alternative hypothesis: two.sided

# Transform data: square root
chloride_primary_sqrt <- Primary_Samples %>%
  mutate(Cl...mmol.L. = sqrt(Cl...mmol.L.))

# Rerun lm model with square root transformed data
chloride_primary_sqrt_lm_int <- lm(Cl...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = chloride_primary_sqrt) 
summary(chloride_primary_sqrt_lm_int) # still not significant

## 
## Call:
## lm(formula = Cl...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, 
##     data = chloride_primary_sqrt)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.298717 -0.081125 -0.007769  0.066405  0.297176 
## 
## Coefficients:
##                                         Estimate Std. Error t value Pr(>|t|)
## (Intercept)                            12.603101   0.035643 353.591   <2e-16
## Pregnant.Or.Atresia.L                  -0.001406   0.050407  -0.028    0.978
## Ambient.Or.OAH.L                        0.042392   0.050407   0.841    0.409
## Pregnant.Or.Atresia.L:Ambient.Or.OAH.L -0.029585   0.071286  -0.415    0.682
##                                           
## (Intercept)                            ***
## Pregnant.Or.Atresia.L                     
## Ambient.Or.OAH.L                          
## Pregnant.Or.Atresia.L:Ambient.Or.OAH.L    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1501 on 24 degrees of freedom
## Multiple R-squared:  0.03267,    Adjusted R-squared:  -0.08825 
## F-statistic: 0.2702 on 3 and 24 DF,  p-value: 0.8462

chloride_primary_sqrt_lm_add <- lm(Cl...mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = chloride_primary_sqrt)
summary(chloride_primary_sqrt_lm_add)

## 
## Call:
## lm(formula = Cl...mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, 
##     data = chloride_primary_sqrt)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.29367 -0.08112 -0.01648  0.07480  0.30222 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           12.59897    0.03365 374.377   <2e-16 ***
## Pregnant.Or.Atresia.L  0.01015    0.04132   0.246    0.808    
## Ambient.Or.OAH.L       0.03432    0.04573   0.751    0.460    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1476 on 25 degrees of freedom
## Multiple R-squared:  0.02573,    Adjusted R-squared:  -0.05221 
## F-statistic: 0.3301 on 2 and 25 DF,  p-value: 0.7219

# Check hows it
chloride_primary_sqrt_lm_res_plot <- simulateResiduals(fittedModel = chloride_primary_sqrt_lm_add)
plot(chloride_primary_sqrt_lm_res_plot) # KS still significant

plotResiduals(chloride_primary_sqrt_lm_add)

testDispersion(chloride_primary_sqrt_lm_add)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.9361, p-value = 0.832
## alternative hypothesis: two.sided

# Residuals of model for shapiro test
chloride_primary_lm_int_res <- residuals(chloride_primary_lm_int) 
chloride_primary_lm_add_res <- residuals(chloride_primary_lm_add)


# Normality test
shapiro.test(chloride_primary_lm_int_res) # normal

## 
##  Shapiro-Wilk normality test
## 
## data:  chloride_primary_lm_int_res
## W = 0.96898, p-value = 0.5534

shapiro.test(chloride_primary_lm_add_res) # normal

## 
##  Shapiro-Wilk normality test
## 
## data:  chloride_primary_lm_add_res
## W = 0.96953, p-value = 0.5681

# Scedasticity test: 
# bartlett's test parametric
bartlett.test(Cl...mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples) 

## 
##  Bartlett test of homogeneity of variances
## 
## data:  Cl...mmol.L. by Pregnant.Or.Atresia
## Bartlett's K-squared = 0.76057, df = 1, p-value = 0.3832

bartlett.test(Cl...mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  Cl...mmol.L. by Ambient.Or.OAH
## Bartlett's K-squared = 3.6573, df = 1, p-value = 0.05582

# Levene's test if nonparametric
leveneTest(Cl...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3  1.2426 0.3162
##       24

leveneTest(Cl...mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.1495 0.7021
##       26

leveneTest(Cl...mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value  Pr(>F)  
## group  1  2.9761 0.09637 .
##       26                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Non-parametric regression test???:

Results: Primary Sample (With both treatments and parturition groups)

• lm interactive model not significant (Pr(>|t|) = 0.678), additive model used

Linear Model: Additive lm param ~ parturition + treatment - No significant relationship detected across groups - Evidence: Multiple R-squared: 0.02517, Adjusted R-squared: -0.05281 F-statistic: 0.3228 on 2 and 25 DF, p-value: 0.7271

• How to know which groups are correlated, and how strongly?
• (test lm individually, next code block)

Warning: DHARMA throws error about quantile deviations (dispersion?).

Correlation: Using sqrt transformed chloride values (additive model) - No significant relationship detected with glucose and parturition groups or treatment groups. - Evidence: Multiple R-squared: 0.02573, Adjusted R-squared: -0.05221 F-statistic: 0.3301 on 2 and 25 DF, p-value: 0.7219

Parametric Assumptions: Across Sample comparisons…

• Sample size (n > 30): Not Met

• Normality (of residuals) : met

• Homoscedasticity: Not Met (for treatment group)

• Evidence:

• Parturition group: Bartlett’s K-squared = 0.76057, df = 1, p-value = 0.3832

• Treatment group: Bartlett’s K-squared = 3.6573, df = 1, p-value = 0.05582

Non-parametric Tests:

• nonparametric normality test (Levene’s test): Met, sorta for treatment group

• Nonparametric regression stat test???

# chloride
# Primary Samples

# Regression Model pt 2

# Individually assess parameter to parturition groups & treatment groups
# p-value significant if (a < 0.1 or a < 0.05)
chloride_primary_lm1 <- lm(Cl...mmol.L. ~ Pregnant.Or.Atresia, data = chloride_primary_sqrt)
summary(chloride_primary_lm1) 

## 
## Call:
## lm(formula = Cl...mmol.L. ~ Pregnant.Or.Atresia, data = chloride_primary_sqrt)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.307152 -0.087238  0.003367  0.068150  0.288741 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           12.58630    0.02887 436.038   <2e-16 ***
## Pregnant.Or.Atresia.L  0.01282    0.04082   0.314    0.756    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1464 on 26 degrees of freedom
## Multiple R-squared:  0.003776,   Adjusted R-squared:  -0.03454 
## F-statistic: 0.09856 on 1 and 26 DF,  p-value: 0.7561

chloride_primary_lm2 <- lm(Cl...mmol.L. ~ Ambient.Or.OAH, data = chloride_primary_sqrt)
summary(chloride_primary_lm2) 

## 
## Call:
## lm(formula = Cl...mmol.L. ~ Ambient.Or.OAH, data = chloride_primary_sqrt)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.28820 -0.08641 -0.01170  0.07993  0.30769 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      12.60136    0.03163 398.457   <2e-16 ***
## Ambient.Or.OAH.L  0.03528    0.04473   0.789    0.437    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1449 on 26 degrees of freedom
## Multiple R-squared:  0.02338,    Adjusted R-squared:  -0.01418 
## F-statistic: 0.6224 on 1 and 26 DF,  p-value: 0.4373

chloride_primary_lm_res_plot1 <- simulateResiduals(fittedModel = chloride_primary_lm1)
plot(chloride_primary_lm_res_plot1)

chloride_primary_lm_res_plot2 <- simulateResiduals(fittedModel = chloride_primary_lm2)
plot(chloride_primary_lm_res_plot2)

visualize(chloride_primary_lm1, plot = "model") 

visualize(chloride_primary_lm2, plot = "model")

# For shapiro test
chloride_primary_lm_res1 <- residuals(chloride_primary_lm1)

chloride_primary_lm_res2 <- residuals(chloride_primary_lm2)

# Normality test
shapiro.test(chloride_primary_lm_res1) # Normal

## 
##  Shapiro-Wilk normality test
## 
## data:  chloride_primary_lm_res1
## W = 0.97367, p-value = 0.6816

shapiro.test(chloride_primary_lm_res2) # Normal

## 
##  Shapiro-Wilk normality test
## 
## data:  chloride_primary_lm_res2
## W = 0.96785, p-value = 0.5243

# Scedasticity test: On sqrt transformed data
bartlett.test(Cl...mmol.L. ~ Pregnant.Or.Atresia, data = chloride_primary_sqrt) 

## 
##  Bartlett test of homogeneity of variances
## 
## data:  Cl...mmol.L. by Pregnant.Or.Atresia
## Bartlett's K-squared = 0.75694, df = 1, p-value = 0.3843

bartlett.test(Cl...mmol.L. ~ Ambient.Or.OAH, data = chloride_primary_sqrt)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  Cl...mmol.L. by Ambient.Or.OAH
## Bartlett's K-squared = 3.6841, df = 1, p-value = 0.05494

# Nonparametric variance test:
leveneTest(Cl...mmol.L. ~ Pregnant.Or.Atresia, data = chloride_primary_sqrt)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.1473 0.7042
##       26

leveneTest(Cl...mmol.L. ~ Ambient.Or.OAH, data = chloride_primary_sqrt)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value  Pr(>F)  
## group  1  3.0146 0.09436 .
##       26                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Results: Primary Sample Regression Analysis pt2

Linear Model: Individual test correlation between param and treatment, and param and parturition outcome

Correlation: Using sqrt transformed sodium values - No significant relationship detected with glucose and parturition groups or treatment groups. - Evidence: - Meaning:

Warning: DHARMA says all good for square root transformed lm model

Parametric Assumptions: Across Sample comparisons…

• Sample size (n > 30): Not Met

• Normality: Normal

• Homoscedasticity: Not Met, for sqrt transformed data, treatment group

• Evidence: On sqrt transformed data

• Parturition group: Bartlett’s K-squared = 0.75694, df = 1, p-value = 0.3843

• Treatment group: Bartlett’s K-squared = 3.6841, df = 1, p-value = 0.05494

• Nonparametric variance test (levene’s test):

• Not met for treatment group: Pr(>F) = 0.09436 .

Residuals Plot (if significant)

# chloride
# Primary samples

# ggplot with geom_smooth(aes(x=,y=)), method = "lm")


# lm for chloride & parturition scatterplot:
chloride_primary_lm_scatterplot1 <- ggplot(data = chloride_primary_sqrt, aes(x = Pregnant.Or.Atresia, y = Cl...mmol.L.)) +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Chloride", x = "Parturition Group", y = "Sqrt Chloride", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(chloride_primary_lm_scatterplot1)

## `geom_smooth()` using formula = 'y ~ x'

ggsave(filename = "data-images/chloride_primary_lm_scatterplot1.pdf", plot = chloride_primary_lm_scatterplot1, device = "pdf")

## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'

ggsave(filename = "data-images/chloride_primary_lm_scatterplot1.png", plot = chloride_primary_lm_scatterplot1, device = "png")

## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'

# lm for glucose & treatment scatterplot:
chloride_primary_lm_scatterplot2 <- ggplot(data = chloride_primary_sqrt, aes(x = Ambient.Or.OAH, y = Cl...mmol.L.)) +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Chloride", x = "Treatment Group", y = "Sqrt Chloride", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(chloride_primary_lm_scatterplot2)

## `geom_smooth()` using formula = 'y ~ x'

ggsave(filename = "data-images/chloride_primary_lm_scatterplot2.pdf", plot = chloride_primary_lm_scatterplot2, device = "pdf")

## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'

ggsave(filename = "data-images/chloride_primary_lm_scatterplot2.png", plot = chloride_primary_lm_scatterplot2, device = "png")

## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'

# chloride
# lm Boxplot

chloride_primary_lm_boxplot1 <- ggplot(data = chloride_primary_sqrt, aes(x = Pregnant.Or.Atresia, y = Cl...mmol.L.)) +
  geom_boxplot() +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Chloride", x = "Parturition Group", y = "Sqrt Chloride", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))
print(chloride_primary_lm_boxplot1)

## `geom_smooth()` using formula = 'y ~ x'

chloride_primary_lm_boxplot2 <- ggplot(data = chloride_primary_sqrt, aes(x = Ambient.Or.OAH, y = Cl...mmol.L.)) +
  geom_boxplot() +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Chloride", x = "Parturition Group", y = "Sqrt Chloride", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))
print(chloride_primary_lm_boxplot2)

## `geom_smooth()` using formula = 'y ~ x'

K+ : Potassium

Potassium Summary Stats

# K+
# Summary stats

# Primary data
Primary_Samples %>%
  group_by(Pregnant.Or.Atresia, Ambient.Or.OAH) %>%
  summarize(count = n(),
            median = round(median(K...mmol.L.), 3),
            mean = round(mean(K...mmol.L.), 3),
            sd = round(sd(K...mmol.L.), 3),
            cv = round(sd(K...mmol.L.)/mean(K...mmol.L.), 3)) %>%
  ungroup()

## `summarise()` has grouped output by 'Pregnant.Or.Atresia'. You can override
## using the `.groups` argument.

## # A tibble: 4 × 7
##   Pregnant.Or.Atresia Ambient.Or.OAH count median  mean    sd    cv
##   <ord>               <ord>          <int>  <dbl> <dbl> <dbl> <dbl>
## 1 Post Parturition    Ambient            8    2.7  2.71 0.146 0.054
## 2 Post Parturition    OAH                2    3    3    0.283 0.094
## 3 Atresia             Ambient           13    2.7  2.86 0.304 0.106
## 4 Atresia             OAH                5    2.7  3    0.714 0.238

Potassium Boxplot

# K+ boxplot

# Primary samples
potassium_primary_boxplot <- ggplot(data = Primary_Samples, aes(x = Pregnant.Or.Atresia, y = K...mmol.L., fill = Ambient.Or.OAH)) +
  geom_boxplot(aes(x = Pregnant.Or.Atresia, y = K...mmol.L., fill = Ambient.Or.OAH)) +
  geom_point(aes(color = Ambient.Or.OAH), position = position_dodge(width = 0.75), alpha = 0.3, color = "black") +
  labs(title = "Potassium", x = "Parturition Type", y = "Potassium (mmol/L)") +
  guides(fill = guide_legend((title = "Treatment Type"))) + 
  scale_fill_manual(values = c("Ambient" = "#189bff", "OAH" = "#598c78")) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(potassium_primary_boxplot)

ggsave(filename = "data-images/potassium_primary_boxplot.pdf", plot = potassium_primary_boxplot, device = "pdf")

## Saving 7 x 5 in image

ggsave(filename = "data-images/potassium_primary_boxplot.png", plot = potassium_primary_boxplot, device = "png")

## Saving 7 x 5 in image

Data Distribution

# K+
# Parametric Assumptions

# Data Distribution
hist(Primary_Samples$K...mmol.L.)

plotNormalHistogram(Primary_Samples$K...mmol.L.) 

plotNormalDensity(Primary_Samples$K...mmol.L.)

K+ Stat Tests

Differences:

# K+
# Primary Samples

# Interactive aov model
potassium_primary_aov_int <- aov(K...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)
summary(potassium_primary_aov_int) # interaction not significant

##                                    Df Sum Sq Mean Sq F value Pr(>F)
## Pregnant.Or.Atresia                 1  0.109 0.10864   0.772  0.388
## Ambient.Or.OAH                      1  0.177 0.17685   1.256  0.274
## Pregnant.Or.Atresia:Ambient.Or.OAH  1  0.025 0.02463   0.175  0.680
## Residuals                          24  3.380 0.14081

# Additive aov model
# If interaction not significant (p-value > 0.05), use additive model
potassium_primary_aov_add <- aov(K...mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = Primary_Samples)
summary(potassium_primary_aov_add)

##                     Df Sum Sq Mean Sq F value Pr(>F)
## Pregnant.Or.Atresia  1  0.109  0.1086   0.798  0.380
## Ambient.Or.OAH       1  0.177  0.1769   1.299  0.265
## Residuals           25  3.404  0.1362

# Parametric Post-Hoc: Tukey's HSD test
TukeyHSD(potassium_primary_aov_add)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = K...mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = Primary_Samples)
## 
## $Pregnant.Or.Atresia
##                          diff        lwr       upr     p adj
## Atresia-Post Parturition 0.13 -0.1697415 0.4297415 0.3802489
## 
## $Ambient.Or.OAH
##                  diff        lwr       upr    p adj
## OAH-Ambient 0.1828571 -0.1488268 0.5145411 0.266967

# Normality plot

# Qqplot: Residuals of aov or lm model (plot residuals of aov vs lm are identical)
potassium_primary_res_qqplot <- ggqqplot(residuals(aov(K...mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = Primary_Samples))) +
labs(title = "Potassium Interactive Residual QQplot",
     subtitle = "residuals(aov(K...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples))",
                               x = "Potassium Theoretical Quantiles (Predicted values)",
                               y = "Potassium Sample Quantiles") +
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(size = 10, hjust = 0.5))
print(potassium_primary_res_qqplot)

# Normality Test
shapiro.test(Primary_Samples$K...mmol.L.)

## 
##  Shapiro-Wilk normality test
## 
## data:  Primary_Samples$K...mmol.L.
## W = 0.76701, p-value = 2.949e-05

shapiro.test(residuals(lm(K...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(lm(K...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples))
## W = 0.86428, p-value = 0.00184

# Parametric variance test:
bartlett.test(K...mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples) 

## 
##  Bartlett test of homogeneity of variances
## 
## data:  K...mmol.L. by Pregnant.Or.Atresia
## Bartlett's K-squared = 5.4475, df = 1, p-value = 0.0196

bartlett.test(K...mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  K...mmol.L. by Ambient.Or.OAH
## Bartlett's K-squared = 7.2218, df = 1, p-value = 0.007202

# Nonparametric variance test:
LeveneTest(K...mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  1.3073 0.2633
##       26

leveneTest(K...mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value  Pr(>F)  
## group  1  3.3217 0.07989 .
##       26                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Nonparametric Stat test & Post-Hoc:

kruskal.test(K...mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  K...mmol.L. by Pregnant.Or.Atresia
## Kruskal-Wallis chi-squared = 0.15497, df = 1, p-value = 0.6938

kruskal.test(K...mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  K...mmol.L. by Ambient.Or.OAH
## Kruskal-Wallis chi-squared = 0.24017, df = 1, p-value = 0.6241

DunnTest(K...mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)

## 
##  Dunn's test of multiple comparisons using rank sums : holm  
## 
##                          mean.rank.diff   pval    
## Atresia-Post Parturition       1.244444 0.6938    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

DunnTest(K...mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)

## 
##  Dunn's test of multiple comparisons using rank sums : holm  
## 
##             mean.rank.diff   pval    
## OAH-Ambient       1.714286 0.6241    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

visualize(potassium_primary_aov_add)

# Test outliers
outlierTest(potassium_primary_aov_add) 

##    rstudent unadjusted p-value Bonferroni p
## 17 4.648451         0.00010152    0.0028426

#  Used Ranked model to adjust for outlier
boxplot(rank(K...mmol.L.) ~ Pregnant.Or.Atresia, data = Primary_Samples, horizontal = TRUE)

boxplot(rank(K...mmol.L.) ~ Ambient.Or.OAH, data = Primary_Samples, horizontal = TRUE)

# Nonparametric Wilcoxon test
wilcox.test(rank(K...mmol.L.) ~ Pregnant.Or.Atresia, data = Primary_Samples)

## Warning in wilcox.test.default(x = DATA[[1L]], y = DATA[[2L]], ...): cannot
## compute exact p-value with ties

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  rank(K...mmol.L.) by Pregnant.Or.Atresia
## W = 82, p-value = 0.7121
## alternative hypothesis: true location shift is not equal to 0

wilcox.test(rank(K...mmol.L.) ~ Ambient.Or.OAH, data = Primary_Samples)

## Warning in wilcox.test.default(x = DATA[[1L]], y = DATA[[2L]], ...): cannot
## compute exact p-value with ties

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  rank(K...mmol.L.) by Ambient.Or.OAH
## W = 64.5, p-value = 0.6435
## alternative hypothesis: true location shift is not equal to 0

Parametric Results: Ambient Fecundity No Atresia Sample

• Interaction not significant (Pr(>F) = 0.680), additive model used.

• Additive summary results:

• No significant differences between parturition and treatment groups.

• Parametric post Hoc: Tukey HSD with additive model

• No groups appear significantly different from eachother…

• Evidence:

Parametric Assumptions: Across Sample comparisons…

• Sample size (n > 30): Not Met

• Normality: Not Met

• Homoscedasticity: Not met (for both factors)

• Evidence:

• Parturition group:

• Treatment group:

Non-Parametric Results:

• nonparametric scedasticity test (Levene’s test): Met ish
• Evidence:
• Parturition group: F value = 1.3073, Pr(>F) = 0.2633
• Treatment group: F value = 3.3217, Pr(>F) = 0.07989

Kruskal-Wallis rank sum test: - No significant difference in chloride detected between parturition and treatment factors - Evidence:

Nonparametric Post-Hoc: Dunn’s test - No significant difference detected across factor groups. - Evidence: - Parturition group: - Treatment group:

Outlier test: Appears significant

Similarities:

# potassium
# Primary Samples

# Regression Model 
# p-value significant if (a < 0.1 or a < 0.05)

potassium_primary_lm_int <- lm(K...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples) # not significant
summary(potassium_primary_lm_int)

## 
## Call:
## lm(formula = K...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, 
##     data = Primary_Samples)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.50000 -0.17115 -0.03702  0.09063  1.20000 
## 
## Coefficients:
##                                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                             2.89351    0.08909  32.477   <2e-16 ***
## Pregnant.Or.Atresia.L                   0.05269    0.12600   0.418    0.680    
## Ambient.Or.OAH.L                        0.15060    0.12600   1.195    0.244    
## Pregnant.Or.Atresia.L:Ambient.Or.OAH.L -0.07452    0.17819  -0.418    0.680    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3753 on 24 degrees of freedom
## Multiple R-squared:  0.08405,    Adjusted R-squared:  -0.03044 
## F-statistic: 0.7341 on 3 and 24 DF,  p-value: 0.5419

potassium_primary_lm_add <- lm(K...mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = Primary_Samples)
summary(potassium_primary_lm_add)

## 
## Call:
## lm(formula = K...mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, 
##     data = Primary_Samples)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.53305 -0.14883 -0.04099  0.06687  1.16695 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            2.88310    0.08412  34.272   <2e-16 ***
## Pregnant.Or.Atresia.L  0.08179    0.10329   0.792    0.436    
## Ambient.Or.OAH.L       0.13026    0.11430   1.140    0.265    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.369 on 25 degrees of freedom
## Multiple R-squared:  0.07738,    Adjusted R-squared:  0.003568 
## F-statistic: 1.048 on 2 and 25 DF,  p-value: 0.3654

# potassium_primary_lm_res_plot <- simulateResiduals(fittedModel = potassium_primary_lm_int)
# plot(potassium_primary_lm_res_plot)

potassium_primary_lm_res_plot <- simulateResiduals(fittedModel = potassium_primary_lm_add)
plot(potassium_primary_lm_res_plot) # quantile deviations detected

visualize(potassium_primary_lm_add, plot = "model")

# Check Dispersion or Outliers
plotResiduals(potassium_primary_lm_add)

testDispersion(potassium_primary_lm_add)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.9361, p-value = 0.832
## alternative hypothesis: two.sided

# For shapiro.test
potassium_primary_lm_add_res <- residuals(potassium_primary_lm_add)

# Normality test
shapiro.test(potassium_primary_lm_add_res) # nonparametric

## 
##  Shapiro-Wilk normality test
## 
## data:  potassium_primary_lm_add_res
## W = 0.87247, p-value = 0.002751

# Scedasticity test: 
# bartlett's test parametric
bartlett.test(K...mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples) 

## 
##  Bartlett test of homogeneity of variances
## 
## data:  K...mmol.L. by Pregnant.Or.Atresia
## Bartlett's K-squared = 5.4475, df = 1, p-value = 0.0196

bartlett.test(K...mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  K...mmol.L. by Ambient.Or.OAH
## Bartlett's K-squared = 7.2218, df = 1, p-value = 0.007202

# Nonparametric variance test (Levene's test):
leveneTest(K...mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  1.3073 0.2633
##       26

leveneTest(K...mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value  Pr(>F)  
## group  1  3.3217 0.07989 .
##       26                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Results: Primary Sample (With both treatments and parturition groups)

• lm interactive model not significant (Pr(>|t|) = 0.680), additive model used

Linear Model: Additive lm param ~ parturition + treatment - No significant relationship detected across groups - Evidence: Multiple R-squared: 0.07738, Adjusted R-squared: 0.003568 F-statistic: 1.048 on 2 and 25 DF, p-value: 0.3654

• How to know which groups are correlated, and how strongly?
• (test lm individually, next code block)

Warning: No DHARMA errors (except a star).

Parametric Assumptions: Across Sample comparisons…

• Sample size (n > 30): Not Met

• Normality (of residuals) : Not met

• Homoscedasticity: Not Met

• Evidence:

• Parturition group:

• Treatment group:

Non-parametric results:

nonparametric variance test (Levene’s test): - Homoscedasticity: Not met for treatment group

transform data (if necessary)

# potassium
# Primary Samples

# Transform data: square root
potassium_primary_sqrt <- Primary_Samples %>%
  mutate(K...mmol.L. = sqrt(K...mmol.L.))

# Rerun lm model with square root transformed data
potassium_primary_sqrt_lm_int <- lm(K...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = potassium_primary_sqrt) 
summary(potassium_primary_sqrt_lm_int) # still not significant

## 
## Call:
## lm(formula = K...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, 
##     data = potassium_primary_sqrt)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.14196 -0.04920 -0.00974  0.02954  0.32629 
## 
## Coefficients:
##                                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                             1.69754    0.02503  67.819   <2e-16 ***
## Pregnant.Or.Atresia.L                   0.01240    0.03540   0.350    0.729    
## Ambient.Or.OAH.L                        0.04180    0.03540   1.181    0.249    
## Pregnant.Or.Atresia.L:Ambient.Or.OAH.L -0.02552    0.05006  -0.510    0.615    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1054 on 24 degrees of freedom
## Multiple R-squared:  0.08134,    Adjusted R-squared:  -0.03349 
## F-statistic: 0.7084 on 3 and 24 DF,  p-value: 0.5564

potassium_primary_sqrt_lm_add <- lm(K...mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = potassium_primary_sqrt)
summary(potassium_primary_sqrt_lm_add)

## 
## Call:
## lm(formula = K...mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, 
##     data = potassium_primary_sqrt)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.15328 -0.04199 -0.01110  0.02141  0.31497 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            1.69398    0.02368  71.549   <2e-16 ***
## Pregnant.Or.Atresia.L  0.02237    0.02907   0.769    0.449    
## Ambient.Or.OAH.L       0.03483    0.03217   1.083    0.289    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1039 on 25 degrees of freedom
## Multiple R-squared:  0.0714, Adjusted R-squared:  -0.002893 
## F-statistic: 0.9611 on 2 and 25 DF,  p-value: 0.3962

# Check hows it
potassium_primary_sqrt_lm_res_plot <- simulateResiduals(fittedModel = potassium_primary_sqrt_lm_add)
plot(potassium_primary_sqrt_lm_res_plot) # KS still significant

plotResiduals(potassium_primary_sqrt_lm_add)

testDispersion(potassium_primary_sqrt_lm_add)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.9361, p-value = 0.832
## alternative hypothesis: two.sided

# Residuals of model for shapiro test
potassium_primary_lm_add_sqrt_res <- residuals(potassium_primary_sqrt_lm_add)

# Normality test
shapiro.test(potassium_primary_lm_add_sqrt_res) # nonparametric

## 
##  Shapiro-Wilk normality test
## 
## data:  potassium_primary_lm_add_sqrt_res
## W = 0.89113, p-value = 0.007134

# Scedasticity test: 
# bartlett's test parametric
bartlett.test(K...mmol.L. ~ Pregnant.Or.Atresia, data = potassium_primary_sqrt) 

## 
##  Bartlett test of homogeneity of variances
## 
## data:  K...mmol.L. by Pregnant.Or.Atresia
## Bartlett's K-squared = 4.737, df = 1, p-value = 0.02952

bartlett.test(K...mmol.L. ~ Ambient.Or.OAH, data = potassium_primary_sqrt)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  K...mmol.L. by Ambient.Or.OAH
## Bartlett's K-squared = 6.3646, df = 1, p-value = 0.01164

# Levene's test if nonparametric
leveneTest(K...mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = potassium_primary_sqrt)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3  1.4427 0.2551
##       24

leveneTest(K...mmol.L. ~ Pregnant.Or.Atresia, data = potassium_primary_sqrt)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  1.3218 0.2607
##       26

leveneTest(K...mmol.L. ~ Ambient.Or.OAH, data = potassium_primary_sqrt)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value  Pr(>F)  
## group  1   3.347 0.07882 .
##       26                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Non-parametric regression test???:

Correlation: Using sqrt transformed chloride values (additive model) - No significant relationship detected with glucose and parturition groups or treatment groups. - Evidence:

Parametric Assumptions: Across Sample comparisons…

• Sample size (n > 30): Not Met

• Normality (of residuals) : Not met

• Homoscedasticity: Not Met (for both groups)

• Evidence:

• Parturition group:

• Treatment group:

Non-parametric Tests:

• nonparametric normality test (Levene’s test): Not really met for treatment group

• Nonparametric regression stat test???

# chloride
# Primary Samples

# Regression Model pt 2

# Individually assess parameter to parturition groups & treatment groups
# p-value significant if (a < 0.1 or a < 0.05)
chloride_primary_lm1 <- lm(Cl...mmol.L. ~ Pregnant.Or.Atresia, data = chloride_primary_sqrt)
summary(chloride_primary_lm1) 

## 
## Call:
## lm(formula = Cl...mmol.L. ~ Pregnant.Or.Atresia, data = chloride_primary_sqrt)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.307152 -0.087238  0.003367  0.068150  0.288741 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           12.58630    0.02887 436.038   <2e-16 ***
## Pregnant.Or.Atresia.L  0.01282    0.04082   0.314    0.756    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1464 on 26 degrees of freedom
## Multiple R-squared:  0.003776,   Adjusted R-squared:  -0.03454 
## F-statistic: 0.09856 on 1 and 26 DF,  p-value: 0.7561

chloride_primary_lm2 <- lm(Cl...mmol.L. ~ Ambient.Or.OAH, data = chloride_primary_sqrt)
summary(chloride_primary_lm2) 

## 
## Call:
## lm(formula = Cl...mmol.L. ~ Ambient.Or.OAH, data = chloride_primary_sqrt)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.28820 -0.08641 -0.01170  0.07993  0.30769 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      12.60136    0.03163 398.457   <2e-16 ***
## Ambient.Or.OAH.L  0.03528    0.04473   0.789    0.437    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1449 on 26 degrees of freedom
## Multiple R-squared:  0.02338,    Adjusted R-squared:  -0.01418 
## F-statistic: 0.6224 on 1 and 26 DF,  p-value: 0.4373

chloride_primary_lm_res_plot1 <- simulateResiduals(fittedModel = chloride_primary_lm1)
plot(chloride_primary_lm_res_plot1)

chloride_primary_lm_res_plot2 <- simulateResiduals(fittedModel = chloride_primary_lm2)
plot(chloride_primary_lm_res_plot2)

visualize(chloride_primary_lm1, plot = "model") 

visualize(chloride_primary_lm2, plot = "model")

# For shapiro test
chloride_primary_lm_res1 <- residuals(chloride_primary_lm1)

chloride_primary_lm_res2 <- residuals(chloride_primary_lm2)

# Normality test
shapiro.test(chloride_primary_lm_res1) # Normal

## 
##  Shapiro-Wilk normality test
## 
## data:  chloride_primary_lm_res1
## W = 0.97367, p-value = 0.6816

shapiro.test(chloride_primary_lm_res2) # Normal

## 
##  Shapiro-Wilk normality test
## 
## data:  chloride_primary_lm_res2
## W = 0.96785, p-value = 0.5243

# Scedasticity test: On sqrt transformed data
bartlett.test(Cl...mmol.L. ~ Pregnant.Or.Atresia, data = chloride_primary_sqrt) 

## 
##  Bartlett test of homogeneity of variances
## 
## data:  Cl...mmol.L. by Pregnant.Or.Atresia
## Bartlett's K-squared = 0.75694, df = 1, p-value = 0.3843

bartlett.test(Cl...mmol.L. ~ Ambient.Or.OAH, data = chloride_primary_sqrt)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  Cl...mmol.L. by Ambient.Or.OAH
## Bartlett's K-squared = 3.6841, df = 1, p-value = 0.05494

# Nonparametric variance test:
leveneTest(Cl...mmol.L. ~ Pregnant.Or.Atresia, data = chloride_primary_sqrt)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.1473 0.7042
##       26

leveneTest(Cl...mmol.L. ~ Ambient.Or.OAH, data = chloride_primary_sqrt)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value  Pr(>F)  
## group  1  3.0146 0.09436 .
##       26                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# potassium
# Primary Samples

# Regression Model pt 2: Not using sqrt model

# Individually assess parameter to parturition groups & treatment groups
# p-value significant if (a < 0.1 or a < 0.05)
potassium_primary_lm1 <- lm(K...mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)
summary(potassium_primary_lm1) 

## 
## Call:
## lm(formula = K...mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.4000 -0.2000 -0.0700  0.0725  1.3000 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            2.83500    0.07319  38.737   <2e-16 ***
## Pregnant.Or.Atresia.L  0.09192    0.10350   0.888    0.383    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3711 on 26 degrees of freedom
## Multiple R-squared:  0.02945,    Adjusted R-squared:  -0.007884 
## F-statistic: 0.7888 on 1 and 26 DF,  p-value: 0.3826

potassium_primary_lm2 <- lm(K...mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)
summary(potassium_primary_lm2) 

## 
## Call:
## lm(formula = K...mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.50000 -0.20119 -0.10476  0.09643  1.20000 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       2.90238    0.07994  36.305   <2e-16 ***
## Ambient.Or.OAH.L  0.13805    0.11306   1.221    0.233    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3664 on 26 degrees of freedom
## Multiple R-squared:  0.05424,    Adjusted R-squared:  0.01786 
## F-statistic: 1.491 on 1 and 26 DF,  p-value: 0.233

potassium_primary_lm_res_plot1 <- simulateResiduals(fittedModel = potassium_primary_lm1)
plot(potassium_primary_lm_res_plot1)

potassium_primary_lm_res_plot2 <- simulateResiduals(fittedModel = potassium_primary_lm2)
plot(potassium_primary_lm_res_plot2)

visualize(potassium_primary_lm1, plot = "model") 

visualize(potassium_primary_lm2, plot = "model")

# For shapiro test
potassium_primary_lm_res1 <- residuals(potassium_primary_lm1)

potassium_primary_lm_res2 <- residuals(potassium_primary_lm2)

# Normality test
shapiro.test(potassium_primary_lm_res1) # Normal

## 
##  Shapiro-Wilk normality test
## 
## data:  potassium_primary_lm_res1
## W = 0.81788, p-value = 0.0002244

shapiro.test(potassium_primary_lm_res2) # Normal

## 
##  Shapiro-Wilk normality test
## 
## data:  potassium_primary_lm_res2
## W = 0.84231, p-value = 0.0006565

# Scedasticity test: On sqrt transformed data
bartlett.test(K...mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples) 

## 
##  Bartlett test of homogeneity of variances
## 
## data:  K...mmol.L. by Pregnant.Or.Atresia
## Bartlett's K-squared = 5.4475, df = 1, p-value = 0.0196

bartlett.test(K...mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  K...mmol.L. by Ambient.Or.OAH
## Bartlett's K-squared = 7.2218, df = 1, p-value = 0.007202

# Nonparametric variance test:
leveneTest(K...mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  1.3073 0.2633
##       26

leveneTest(K...mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value  Pr(>F)  
## group  1  3.3217 0.07989 .
##       26                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residuals Plot (if significant)

# potassium
# Primary samples

# ggplot with geom_smooth(aes(x=,y=)), method = "lm")


# lm for chloride & parturition scatterplot:
potassium_primary_lm_scatterplot1 <- ggplot(data = Primary_Samples, aes(x = Pregnant.Or.Atresia, y = K...mmol.L.)) +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Potassium", x = "Parturition Group", y = "Potassium (mmol/L)", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(potassium_primary_lm_scatterplot1)

## `geom_smooth()` using formula = 'y ~ x'

ggsave(filename = "data-images/potassium_primary_lm_scatterplot1.pdf", plot = potassium_primary_lm_scatterplot1, device = "pdf")

## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'

ggsave(filename = "data-images/potassium_primary_lm_scatterplot1.png", plot = potassium_primary_lm_scatterplot1, device = "png")

## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'

# lm for glucose & treatment scatterplot:
potassium_primary_lm_scatterplot2 <- ggplot(data = Primary_Samples, aes(x = Ambient.Or.OAH, y = K...mmol.L.)) +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Potassium", x = "Treatment Group", y = "Potassium (mmol/L)", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(potassium_primary_lm_scatterplot2)

## `geom_smooth()` using formula = 'y ~ x'

ggsave(filename = "data-images/potassium_primary_lm_scatterplot2.pdf", plot = potassium_primary_lm_scatterplot2, device = "pdf")

## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'

ggsave(filename = "data-images/potassium_primary_lm_scatterplot2.png", plot = potassium_primary_lm_scatterplot2, device = "png")

## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'

# K+
# lm Boxplot

potassium_primary_lm_boxplot1 <- ggplot(data = Primary_Samples, aes(x = Pregnant.Or.Atresia, y = K...mmol.L.)) +
  geom_boxplot() +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Potassium", x = "Parturition Group", y = "Potassium (mmol/L)", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))
print(potassium_primary_lm_boxplot1)

## `geom_smooth()` using formula = 'y ~ x'

potassium_primary_lm_boxplot2 <- ggplot(data = Primary_Samples, aes(x = Ambient.Or.OAH, y = K...mmol.L.)) +
  geom_boxplot() +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Potassium", x = "Parturition Group", y = "Potassium (mmol/L)", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))
print(potassium_primary_lm_boxplot2)

## `geom_smooth()` using formula = 'y ~ x'

Ca+2 : Calcium

Calcium Summary Stats

# Ca+2
# Summary stats

# Primary data
Primary_Samples %>%
  group_by(Pregnant.Or.Atresia, Ambient.Or.OAH) %>%
  summarize(count = n(),
            median = round(median(Ca....mmol.L.), 3),
            mean = round(mean(Ca....mmol.L.), 3),
            sd = round(sd(Ca....mmol.L.), 3),
            cv = round(sd(Ca....mmol.L.)/mean(Ca....mmol.L.), 3)) %>%
  ungroup()

## `summarise()` has grouped output by 'Pregnant.Or.Atresia'. You can override
## using the `.groups` argument.

## # A tibble: 4 × 7
##   Pregnant.Or.Atresia Ambient.Or.OAH count median  mean    sd    cv
##   <ord>               <ord>          <int>  <dbl> <dbl> <dbl> <dbl>
## 1 Post Parturition    Ambient            8   1.43  1.43 0.041 0.029
## 2 Post Parturition    OAH                2   1.25  1.25 0.071 0.057
## 3 Atresia             Ambient           13   1.24  1.27 0.114 0.09 
## 4 Atresia             OAH                5   1.16  1.16 0.031 0.027

Calcium Boxplot

# Ca+2 boxplot

# Primary samples
calcium_primary_boxplot <- ggplot(data = Primary_Samples, aes(x = Pregnant.Or.Atresia, y = Ca....mmol.L., fill = Ambient.Or.OAH)) +
  geom_boxplot(aes(x = Pregnant.Or.Atresia, y = Ca....mmol.L., fill = Ambient.Or.OAH)) +
  geom_point(aes(color = Ambient.Or.OAH), position = position_dodge(width = 0.75), alpha = 0.3, color = "black") +
  labs(title = "Calcium", x = "Parturition Type", y = "Calcium (mmol/L)") +
  guides(fill = guide_legend((title = "Treatment Type"))) + 
  scale_fill_manual(values = c("Ambient" = "#189bff", "OAH" = "#598c78")) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(calcium_primary_boxplot)

ggsave(filename = "data-images/calcium_primary_boxplot.pdf", plot = calcium_primary_boxplot, device = "pdf")

## Saving 7 x 5 in image

ggsave(filename = "data-images/calcium_primary_boxplot.png", plot = calcium_primary_boxplot, device = "png")

## Saving 7 x 5 in image

Data Distribution

# Ca+2
# Parametric Assumptions

# Data Distribution
hist(Primary_Samples$Ca....mmol.L.)

plotNormalHistogram(Primary_Samples$Ca....mmol.L.) 

plotNormalDensity(Primary_Samples$Ca....mmol.L.)

Ca+2 Stat Tests

Differences:

# Ca+2
# Primary Samples

# Interactive aov model
calcium_primary_aov_int <- aov(Ca....mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)
summary(calcium_primary_aov_int) # interaction not significant

##                                    Df  Sum Sq Mean Sq F value   Pr(>F)    
## Pregnant.Or.Atresia                 1 0.15667 0.15667  21.391 0.000108 ***
## Ambient.Or.OAH                      1 0.08990 0.08990  12.274 0.001827 ** 
## Pregnant.Or.Atresia:Ambient.Or.OAH  1 0.00575 0.00575   0.785 0.384354    
## Residuals                          24 0.17578 0.00732                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Additive aov model
# If interaction not significant (p-value > 0.05), use additive model
calcium_primary_aov_add <- aov(Ca....mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = Primary_Samples)
summary(calcium_primary_aov_add)

##                     Df Sum Sq Mean Sq F value   Pr(>F)    
## Pregnant.Or.Atresia  1 0.1567 0.15667   21.58 9.35e-05 ***
## Ambient.Or.OAH       1 0.0899 0.08990   12.38  0.00169 ** 
## Residuals           25 0.1815 0.00726                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Parametric Post-Hoc: Tukey's HSD test
TukeyHSD(calcium_primary_aov_add)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Ca....mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = Primary_Samples)
## 
## $Pregnant.Or.Atresia
##                                diff        lwr         upr    p adj
## Atresia-Post Parturition -0.1561111 -0.2253289 -0.08689334 9.35e-05
## 
## $Ambient.Or.OAH
##                   diff        lwr         upr     p adj
## OAH-Ambient -0.1303704 -0.2069644 -0.05377629 0.0017412

# Normality plot

# Qqplot: Residuals of aov or lm model (plot residuals of aov vs lm are identical)
calcium_primary_res_qqplot <- ggqqplot(residuals(aov(Ca....mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = Primary_Samples))) +
labs(title = "Calcium Interactive Residual QQplot",
     subtitle = "residuals(aov(Ca....mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples))",
                               x = "Calcium Theoretical Quantiles (Predicted values)",
                               y = "Calcium Sample Quantiles") +
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(size = 10, hjust = 0.5))
print(calcium_primary_res_qqplot)

# Normality Test
shapiro.test(Primary_Samples$Ca....mmol.L.)

## 
##  Shapiro-Wilk normality test
## 
## data:  Primary_Samples$Ca....mmol.L.
## W = 0.93132, p-value = 0.06655

shapiro.test(residuals(lm(Ca....mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples)))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(lm(Ca....mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples))
## W = 0.91748, p-value = 0.03007

# Parametric variance test:
bartlett.test(Ca....mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples) 

## 
##  Bartlett test of homogeneity of variances
## 
## data:  Ca....mmol.L. by Pregnant.Or.Atresia
## Bartlett's K-squared = 0.49921, df = 1, p-value = 0.4798

bartlett.test(Ca....mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  Ca....mmol.L. by Ambient.Or.OAH
## Bartlett's K-squared = 3.5462, df = 1, p-value = 0.05968

# Nonparametric variance test:
LeveneTest(Ca....mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.2268 0.6379
##       26

leveneTest(Ca....mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)  
## group  1   6.801 0.0149 *
##       26                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Nonparametric Stat test & Post-Hoc:

kruskal.test(Ca....mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  Ca....mmol.L. by Pregnant.Or.Atresia
## Kruskal-Wallis chi-squared = 9.7286, df = 1, p-value = 0.001814

kruskal.test(Ca....mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  Ca....mmol.L. by Ambient.Or.OAH
## Kruskal-Wallis chi-squared = 8.0703, df = 1, p-value = 0.0045

DunnTest(Ca....mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)

## 
##  Dunn's test of multiple comparisons using rank sums : holm  
## 
##                          mean.rank.diff   pval    
## Atresia-Post Parturition      -10.11111 0.0018 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

DunnTest(Ca....mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)

## 
##  Dunn's test of multiple comparisons using rank sums : holm  
## 
##             mean.rank.diff   pval    
## OAH-Ambient      -10.19048 0.0045 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

boxplot(Ca....mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples, horizontal = TRUE)

boxplot(Ca....mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples, horizontal = TRUE)

# Test outliers
outlierTest(calcium_primary_aov_add) 

##    rstudent unadjusted p-value Bonferroni p
## 16 3.869496         0.00073227     0.020504

#  Used Ranked model to adjust for outlier
boxplot(rank(Ca....mmol.L.) ~ Pregnant.Or.Atresia, data = Primary_Samples, horizontal = TRUE)

boxplot(rank(Ca....mmol.L.) ~ Ambient.Or.OAH, data = Primary_Samples, horizontal = TRUE)

# Nonparametric Wilcoxon test
wilcox.test(rank(Ca....mmol.L.) ~ Pregnant.Or.Atresia, data = Primary_Samples, alternative = c("two.sided"))

## Warning in wilcox.test.default(x = DATA[[1L]], y = DATA[[2L]], ...): cannot
## compute exact p-value with ties

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  rank(Ca....mmol.L.) by Pregnant.Or.Atresia
## W = 155, p-value = 0.001968
## alternative hypothesis: true location shift is not equal to 0

wilcox.test(rank(Ca....mmol.L.) ~ Ambient.Or.OAH, data = Primary_Samples, alternative = c("two.sided"))

## Warning in wilcox.test.default(x = DATA[[1L]], y = DATA[[2L]], ...): cannot
## compute exact p-value with ties

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  rank(Ca....mmol.L.) by Ambient.Or.OAH
## W = 127, p-value = 0.004889
## alternative hypothesis: true location shift is not equal to 0

Interaction not significant, additive model used.

Results: summary aov additive model - Post parturition group significantly different from atresia group (Pr(>F) = 9.35e-05 *) - Ambient treatment significantly different from OAH treatment (Pr(>F) = 0.00169 )

Parametric Post Hoc: Tukey HSD test - Post parturition patients display a significant increase in calcium levels compared to atresia patients (p adj = 9.35e-05) - Ambient treatment patients display a significant increase in calcium levels compared to OAH treatment patients (p adj = 0.0017412)

Parametric Assumptions:

Normality: Not met - W = 0.93132, p-value = 0.06655

Parametric Scedasticity test (Bartlett’s test): - Not met for treatment groups

Nonparametric variance test: Levene’s test: - Not met for treatment group

Nonparametric Stat test results: Kruskall Wallis test - Post parturition patients display a significant increase in calcium levels compared to atresia patients (Kruskal-Wallis chi-squared = 9.7286, df = 1, p-value = 0.001814) - Ambient treatment patients display a significant increase in calcium levels compared to OAH treatment patients (Kruskal-Wallis chi-squared = 8.0703, df = 1, p-value = 0.0045)

Nonparametric Post Hoc: Dunn’s Test - nparametric Stat test results: Kruskall Wallis test - Post parturition patients display a significant increase in calcium levels compared to atresia patients (OAH-Ambient: mean.rank.diff = -10.19048, pval = 0.0018 ) - Ambient treatment patients display a significant increase in calcium levels compared to OAH treatment patients (OAH-Ambient: mean.rank.diff = -10.19048, pval = 0.0045 )

Similarities:

# calcium
# Primary Samples

# Regression Model 
# p-value significant if (a < 0.1 or a < 0.05)

calcium_primary_lm_int <- lm(Ca....mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = Primary_Samples) # not significant
summary(calcium_primary_lm_int)

## 
## Call:
## lm(formula = Ca....mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, 
##     data = Primary_Samples)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.14923 -0.05000 -0.00524  0.03219  0.26077 
## 
## Coefficients:
##                                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                             1.27762    0.02032  62.878  < 2e-16 ***
## Pregnant.Or.Atresia.L                  -0.08910    0.02874  -3.101  0.00488 ** 
## Ambient.Or.OAH.L                       -0.10270    0.02874  -3.574  0.00153 ** 
## Pregnant.Or.Atresia.L:Ambient.Or.OAH.L  0.03601    0.04064   0.886  0.38435    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.08558 on 24 degrees of freedom
## Multiple R-squared:  0.5894, Adjusted R-squared:  0.5381 
## F-statistic: 11.48 on 3 and 24 DF,  p-value: 7.291e-05

calcium_primary_lm_add <- lm(Ca....mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = Primary_Samples)
summary(calcium_primary_lm_add)

## 
## Call:
## lm(formula = Ca....mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, 
##     data = Primary_Samples)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.155373 -0.042295 -0.003321  0.038470  0.254627 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            1.28265    0.01943  66.026  < 2e-16 ***
## Pregnant.Or.Atresia.L -0.10316    0.02385  -4.325 0.000214 ***
## Ambient.Or.OAH.L      -0.09287    0.02640  -3.519 0.001685 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.08521 on 25 degrees of freedom
## Multiple R-squared:  0.576,  Adjusted R-squared:  0.542 
## F-statistic: 16.98 on 2 and 25 DF,  p-value: 2.201e-05

# calcium_primary_lm_res_plot <- simulateResiduals(fittedModel = calcium_primary_lm_int)
# plot(calcium_primary_lm_res_plot)

calcium_primary_lm_res_plot <- simulateResiduals(fittedModel = calcium_primary_lm_add)
plot(calcium_primary_lm_res_plot) # quantile deviations detected

visualize(calcium_primary_lm_add, plot = "model")

# Check Dispersion or Outliers
plotResiduals(calcium_primary_lm_add)

testDispersion(calcium_primary_lm_add)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.9361, p-value = 0.832
## alternative hypothesis: two.sided

# For shapiro.test
calcium_primary_lm_add_res <- residuals(calcium_primary_lm_add)

# Normality test
shapiro.test(calcium_primary_lm_add_res) # nonparametric

## 
##  Shapiro-Wilk normality test
## 
## data:  calcium_primary_lm_add_res
## W = 0.94439, p-value = 0.1429

# Scedasticity test: 
# bartlett's test parametric
bartlett.test(Ca....mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples) 

## 
##  Bartlett test of homogeneity of variances
## 
## data:  Ca....mmol.L. by Pregnant.Or.Atresia
## Bartlett's K-squared = 0.49921, df = 1, p-value = 0.4798

bartlett.test(Ca....mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  Ca....mmol.L. by Ambient.Or.OAH
## Bartlett's K-squared = 3.5462, df = 1, p-value = 0.05968

# Nonparametric variance test (Levene's test):
leveneTest(Ca....mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.2268 0.6379
##       26

leveneTest(Ca....mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)  
## group  1   6.801 0.0149 *
##       26                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interaction not significant, additive model used.

lm summary results: additive model - Significant relationship detected. - Multiple R-squared: 0.576, Adjusted R-squared: 0.542 F-statistic: 16.98 on 2 and 25 DF, p-value: 2.201e-05

• part 2: for determining which groups have strong correlation

Parametric Assumptions:

Normality: Not met

Parametric Scedasticity test (Bartlett’s test): - Not met for treatment groups

Nonparametric variance test: Levene’s test: - Not met for treatment group

# Ca+2
# Primary Samples

# Regression Model pt 2: Not using sqrt model

# Individually assess parameter to parturition groups & treatment groups
# p-value significant if (a < 0.1 or a < 0.05)
calcium_primary_lm1 <- lm(Ca....mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)
summary(calcium_primary_lm1) 

## 
## Call:
## lm(formula = Ca....mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.195000 -0.063889 -0.006944  0.040000  0.291111 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            1.31694    0.02015  65.360  < 2e-16 ***
## Pregnant.Or.Atresia.L -0.11039    0.02849  -3.874 0.000649 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1022 on 26 degrees of freedom
## Multiple R-squared:  0.366,  Adjusted R-squared:  0.3416 
## F-statistic: 15.01 on 1 and 26 DF,  p-value: 0.000649

calcium_primary_lm2 <- lm(Ca....mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)
summary(calcium_primary_lm2) 

## 
## Call:
## lm(formula = Ca....mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.210952 -0.083452  0.001667  0.091548  0.199048 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       1.25833    0.02411  52.194  < 2e-16 ***
## Ambient.Or.OAH.L -0.10270    0.03409  -3.012  0.00571 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1105 on 26 degrees of freedom
## Multiple R-squared:  0.2587, Adjusted R-squared:  0.2302 
## F-statistic: 9.073 on 1 and 26 DF,  p-value: 0.005714

calcium_primary_lm_res_plot1 <- simulateResiduals(fittedModel = calcium_primary_lm1)
plot(calcium_primary_lm_res_plot1)

calcium_primary_lm_res_plot2 <- simulateResiduals(fittedModel = calcium_primary_lm2)
plot(calcium_primary_lm_res_plot2)

visualize(calcium_primary_lm1, plot = "model") 

visualize(calcium_primary_lm2, plot = "model")

# Check Dispersion or Outliers
plotResiduals(calcium_primary_lm1)

testDispersion(calcium_primary_lm1)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.97355, p-value = 0.928
## alternative hypothesis: two.sided

plotResiduals(calcium_primary_lm2)

testDispersion(calcium_primary_lm2)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.97355, p-value = 0.928
## alternative hypothesis: two.sided

# For shapiro test
calcium_primary_lm_res1 <- residuals(potassium_primary_lm1)

calcium_primary_lm_res2 <- residuals(potassium_primary_lm2)

# Normality test
shapiro.test(calcium_primary_lm_res1) # Normal

## 
##  Shapiro-Wilk normality test
## 
## data:  calcium_primary_lm_res1
## W = 0.81788, p-value = 0.0002244

shapiro.test(calcium_primary_lm_res2) # Normal

## 
##  Shapiro-Wilk normality test
## 
## data:  calcium_primary_lm_res2
## W = 0.84231, p-value = 0.0006565

# Scedasticity test: On sqrt transformed data
bartlett.test(Ca....mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples) 

## 
##  Bartlett test of homogeneity of variances
## 
## data:  Ca....mmol.L. by Pregnant.Or.Atresia
## Bartlett's K-squared = 0.49921, df = 1, p-value = 0.4798

bartlett.test(Ca....mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  Ca....mmol.L. by Ambient.Or.OAH
## Bartlett's K-squared = 3.5462, df = 1, p-value = 0.05968

# Nonparametric variance test:
leveneTest(Ca....mmol.L. ~ Pregnant.Or.Atresia, data = Primary_Samples)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.2268 0.6379
##       26

leveneTest(Ca....mmol.L. ~ Ambient.Or.OAH, data = Primary_Samples)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)  
## group  1   6.801 0.0149 *
##       26                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Results: - A significant correlation detected, increasing calcium levels associated with post-parturition patients opposed to atresia patients.

• A significant correlation detected, increasing calcium levels assoiciated with ambient treatment patients opposed to OAH treatment patients.

Warning DHAMA error thrown for treatment group: Levene test for homogeneity of variance significant

Normality: Not met (for both factor groups)

Parametric Equal Variance test: Bartlett’s test - Not met for treatment group

Nonparametric variance test: Levene’s test - Not met for treatment group

Residuals Plot (if significant)

# Ca+2
# Primary samples

# ggplot with geom_smooth(aes(x=,y=)), method = "lm")


# lm for chloride & parturition scatterplot:
calcium_primary_lm_scatterplot1 <- ggplot(data = Primary_Samples, aes(x = Pregnant.Or.Atresia, y = Ca....mmol.L.)) +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Calcium", x = "Parturition Group", y = "Calcium (mmol/L)", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(calcium_primary_lm_scatterplot1)

## `geom_smooth()` using formula = 'y ~ x'

ggsave(filename = "data-images/calcium_primary_lm_scatterplot1.pdf", plot = calcium_primary_lm_scatterplot1, device = "pdf")

## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'

ggsave(filename = "data-images/calcium_primary_lm_scatterplot1.png", plot = calcium_primary_lm_scatterplot1, device = "png")

## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'

# lm for glucose & treatment scatterplot:
calcium_primary_lm_scatterplot2 <- ggplot(data = Primary_Samples, aes(x = Ambient.Or.OAH, y = Ca....mmol.L.)) +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Calcium", x = "Treatment Group", y = "Calcium (mmol/L)", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(calcium_primary_lm_scatterplot2)

## `geom_smooth()` using formula = 'y ~ x'

ggsave(filename = "data-images/calcium_primary_lm_scatterplot2.pdf", plot = calcium_primary_lm_scatterplot2, device = "pdf")

## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'

ggsave(filename = "data-images/calcium_primary_lm_scatterplot2.png", plot = calcium_primary_lm_scatterplot2, device = "png")

## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'

# Ca+2
# lm Boxplot

calcium_primary_lm_boxplot1 <- ggplot(data = Primary_Samples, aes(x = Pregnant.Or.Atresia, y = Ca....mmol.L.)) +
  geom_boxplot() +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Calcium", x = "Parturition Group", y = "Calcium (mmol/L)", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))
print(calcium_primary_lm_boxplot1)

## `geom_smooth()` using formula = 'y ~ x'

calcium_primary_lm_boxplot2 <- ggplot(data = Primary_Samples, aes(x = Ambient.Or.OAH, y = Ca....mmol.L.)) +
  geom_boxplot() +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Calcium", x = "Parturition Group", y = "Calcium (mmol/L)", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))
print(calcium_primary_lm_boxplot2)

## `geom_smooth()` using formula = 'y ~ x'

transform data (if necessary)

# Ca+2
# Primary Samples

# Transform data: square root
calcium_primary_sqrt <- Primary_Samples %>%
  mutate(Ca....mmol.L. = sqrt(Ca....mmol.L.))

# Rerun lm model with square root transformed data
calcium_primary_sqrt_lm_int <- lm(Ca....mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = calcium_primary_sqrt) 
summary(calcium_primary_sqrt_lm_int) # still not significant

## 
## Call:
## lm(formula = Ca....mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, 
##     data = calcium_primary_sqrt)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.067294 -0.021721 -0.001756  0.014485  0.111337 
## 
## Coefficients:
##                                         Estimate Std. Error t value Pr(>|t|)
## (Intercept)                             1.129150   0.008855 127.513  < 2e-16
## Pregnant.Or.Atresia.L                  -0.039421   0.012523  -3.148  0.00436
## Ambient.Or.OAH.L                       -0.044925   0.012523  -3.587  0.00148
## Pregnant.Or.Atresia.L:Ambient.Or.OAH.L  0.014896   0.017710   0.841  0.40860
##                                           
## (Intercept)                            ***
## Pregnant.Or.Atresia.L                  ** 
## Ambient.Or.OAH.L                       ** 
## Pregnant.Or.Atresia.L:Ambient.Or.OAH.L    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0373 on 24 degrees of freedom
## Multiple R-squared:  0.5924, Adjusted R-squared:  0.5415 
## F-statistic: 11.63 on 3 and 24 DF,  p-value: 6.689e-05

calcium_primary_sqrt_lm_add <- lm(Ca....mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, data = calcium_primary_sqrt)
summary(calcium_primary_sqrt_lm_add)

## 
## Call:
## lm(formula = Ca....mmol.L. ~ Pregnant.Or.Atresia + Ambient.Or.OAH, 
##     data = calcium_primary_sqrt)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.069835 -0.017804 -0.000962  0.017027  0.108797 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            1.131231   0.008453 133.829  < 2e-16 ***
## Pregnant.Or.Atresia.L -0.045238   0.010379  -4.359 0.000196 ***
## Ambient.Or.OAH.L      -0.040860   0.011485  -3.558 0.001527 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03708 on 25 degrees of freedom
## Multiple R-squared:  0.5804, Adjusted R-squared:  0.5468 
## F-statistic: 17.29 on 2 and 25 DF,  p-value: 1.93e-05

# Check hows it
calcium_primary_sqrt_lm_res_plot <- simulateResiduals(fittedModel = calcium_primary_sqrt_lm_add)
plot(calcium_primary_sqrt_lm_res_plot) # KS still significant

plotResiduals(calcium_primary_sqrt_lm_add)

testDispersion(calcium_primary_sqrt_lm_add)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.9361, p-value = 0.832
## alternative hypothesis: two.sided

# Residuals of model for shapiro test
calcium_primary_lm_add_sqrt_res <- residuals(calcium_primary_sqrt_lm_add)

# Normality test
shapiro.test(calcium_primary_lm_add_sqrt_res) # nonparametric

## 
##  Shapiro-Wilk normality test
## 
## data:  calcium_primary_lm_add_sqrt_res
## W = 0.94977, p-value = 0.1955

# Scedasticity test: 
# bartlett's test parametric
bartlett.test(Ca....mmol.L. ~ Pregnant.Or.Atresia, data = calcium_primary_sqrt) 

## 
##  Bartlett test of homogeneity of variances
## 
## data:  Ca....mmol.L. by Pregnant.Or.Atresia
## Bartlett's K-squared = 0.57932, df = 1, p-value = 0.4466

bartlett.test(Ca....mmol.L. ~ Ambient.Or.OAH, data = calcium_primary_sqrt)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  Ca....mmol.L. by Ambient.Or.OAH
## Bartlett's K-squared = 3.1949, df = 1, p-value = 0.07387

# Levene's test if nonparametric
leveneTest(Ca....mmol.L. ~ Pregnant.Or.Atresia * Ambient.Or.OAH, data = calcium_primary_sqrt)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3  1.7959 0.1749
##       24

leveneTest(Ca....mmol.L. ~ Pregnant.Or.Atresia, data = calcium_primary_sqrt)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1   0.327 0.5723
##       26

leveneTest(Ca....mmol.L. ~ Ambient.Or.OAH, data = calcium_primary_sqrt)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value  Pr(>F)  
## group  1  6.1011 0.02039 *
##       26                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Non-parametric regression test???: